mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Compare commits

base: mulimoen:11807ffb684764a305f25b2f579c8603b4bc85e3
Branches Tags
mulimoen:master
..
compare: mulimoen:cfeb30fac02b7e8778e68ee44215ef1152c5571e
Branches Tags
mulimoen:master

No commits in common. "11807ffb684764a305f25b2f579c8603b4bc85e3" and "cfeb30fac02b7e8778e68ee44215ef1152c5571e" have entirely different histories.
11807ffb68 ... cfeb30fac0

37 changed files with 168 additions and 188 deletions
Show Stats Expand all files Collapse all files

12
euler/Cargo.toml
View File

@ -2,7 +2,7 @@
name = "euler"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[features]
# Internal feature flag to gate the expensive tests
@ -11,15 +11,15 @@ expensive_tests = []
serde1 = ["serde", "arrayvec/serde"]
[dependencies]
ndarray = "0.15.4"
ndarray = "0.15.0"
sbp = { path = "../sbp" }
arrayvec = "0.7.2"
serde = { version = "1.0.138", default-features = false, optional = true, features = ["derive"] }
arrayvec = "0.6.0"
serde = { version = "1.0.115", default-features = false, optional = true, features = ["derive"] }
integrate = { path = "../utils/integrate" }
once_cell = "1.13.0"
once_cell = "1.7.2"
[dev-dependencies]
criterion = "0.4.0"
criterion = "0.3.2"
[[bench]]
name = "bench"

48
euler/src/lib.rs
View File

@ -60,7 +60,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
let metrics = &self.grid.1;
let rhs_trad = |k: &mut Diff, y: &Field, _time: Float| {
let boundaries = boundary_extractor(y, grid, &bc);
RHS_trad(op, k, y, metrics, &boundaries, wb);
RHS_trad(op, k, y, metrics, &boundaries, wb)
};
integrate::integrate::<integrate::Rk4, Field, _>(
rhs_trad,
@ -98,7 +98,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
self.sys
.0
.vortex(self.grid.0.x(), self.grid.0.y(), 0.0, &vortex_parameters);
.vortex(self.grid.0.x(), self.grid.0.y(), 0.0, &vortex_parameters)
}
pub fn field(&self) -> &Field {
@ -134,7 +134,7 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
let rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
let (grid, metrics) = grid;
let boundaries = boundary_extractor(y, grid, &bc);
RHS_upwind(op, k, y, metrics, &boundaries, wb);
RHS_upwind(op, k, y, metrics, &boundaries, wb)
};
integrate::integrate::<integrate::Rk4, Field, _>(
rhs_upwind,
@ -159,7 +159,7 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
let mut rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
let (grid, metrics) = grid;
let boundaries = boundary_extractor(y, grid, &bc);
RHS_upwind(op, k, y, metrics, &boundaries, wb);
RHS_upwind(op, k, y, metrics, &boundaries, wb)
};
let mut time = 0.0;
let mut sys2 = self.sys.0.clone();
@ -194,7 +194,7 @@ impl Clone for Field {
Self(self.0.clone())
}
fn clone_from(&mut self, source: &Self) {
self.0.clone_from(&source.0);
self.0.clone_from(&source.0)
}
}
@ -324,7 +324,7 @@ impl Field {
vortex_param: &VortexParameters,
) {
let (rho, rhou, rhov, e) = self.components_mut();
vortex_param.evaluate(time, x, y, rho, rhou, rhov, e);
vortex_param.evaluate(time, x, y, rho, rhou, rhov, e)
}
#[allow(clippy::erasing_op, clippy::identity_op)]
fn iter(&self) -> impl ExactSizeIterator<Item = FieldValue> + '_ {
@ -508,7 +508,7 @@ pub fn RHS_trad(
eflux in &dE.0,
fflux in &dF.0,
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
*out = (-eflux - fflux)/detj;
*out = (-eflux - fflux)/detj
});
SAT_characteristics(op, k, y, metrics, boundaries);
@ -555,14 +555,14 @@ pub fn RHS_no_SAT(
ad_xi in &ad_xi.0,
ad_eta in &ad_eta.0,
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
*out = (-eflux - fflux + ad_xi + ad_eta)/detj;
*out = (-eflux - fflux + ad_xi + ad_eta)/detj
});
} else {
azip!((out in &mut k.0,
eflux in &dE.0,
fflux in &dF.0,
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
*out = (-eflux - fflux )/detj;
*out = (-eflux - fflux )/detj
});
}
}
@ -609,7 +609,7 @@ pub fn RHS_upwind(
ad_xi in &ad_xi.0,
ad_eta in &ad_eta.0,
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
*out = (-eflux - fflux + ad_xi + ad_eta)/detj;
*out = (-eflux - fflux + ad_xi + ad_eta)/detj
});
SAT_characteristics(op, k, y, metrics, boundaries);
@ -687,7 +687,7 @@ fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics, wb: &mut Fi
let mut p = wb.rho_mut();
azip!((p in &mut p, &rho in &rho, &rhou in &rhou, &rhov in &rhov, &e in &e) {
*p = pressure(gamma, rho, rhou, rhov, e);
*p = pressure(gamma, rho, rhou, rhov, e)
});
let (mut c0, c1, mut c2, c3) = k.0.components_mut();
@ -730,7 +730,7 @@ fn fluxes(k: (&mut Field, &mut Field), y: &Field, metrics: &Metrics, wb: &mut Fi
let fflux = *ff;
*ef = j_dxi_dx * eflux + j_dxi_dy * fflux;
*ff = j_deta_dx * eflux + j_deta_dy * fflux;
});
})
}
}
@ -748,11 +748,11 @@ impl std::fmt::Debug for BoundaryCharacteristic {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Self::This => write!(f, "This"),
Self::Grid(g) => write!(f, "Grid({g})"),
Self::Vortex(vp) => write!(f, "{vp:?}"),
Self::Grid(g) => write!(f, "Grid({})", g),
Self::Vortex(vp) => write!(f, "{:?}", vp),
Self::Eval(_) => write!(f, "Eval"),
Self::Interpolate(_, _) => write!(f, "Interpolate"),
Self::MultiGrid(m) => write!(f, "Multigrid: {m:?}"),
Self::MultiGrid(m) => write!(f, "Multigrid: {:?}", m),
}
}
}
@ -835,12 +835,12 @@ fn boundary_extract<'a>(
to.slice_mut(s![.., i..i + n])
.assign(&seldir(&fields[g]).slice(s![.., start..end]));
remaining -= 1;
if remaining == 0 {
i += n;
assert_eq!(i, to.len_of(Axis(1)));
} else {
if remaining != 0 {
to.slice_mut(s![.., i]).iter_mut().for_each(|x| *x /= 2.0);
i += n - 1;
} else {
i += n;
assert_eq!(i, to.len_of(Axis(1)));
}
}
to.view()
@ -880,7 +880,7 @@ pub fn boundary_extracts<'a>(
&bt.north,
field,
grid.north(),
Field::south,
|f| f.south(),
eb.north.as_mut(),
time,
),
@ -889,7 +889,7 @@ pub fn boundary_extracts<'a>(
&bt.south,
field,
grid.south(),
Field::north,
|f| f.north(),
eb.south.as_mut(),
time,
),
@ -898,7 +898,7 @@ pub fn boundary_extracts<'a>(
&bt.east,
field,
grid.east(),
Field::west,
|f| f.west(),
eb.east.as_mut(),
time,
),
@ -907,7 +907,7 @@ pub fn boundary_extracts<'a>(
&bt.west,
field,
grid.west(),
Field::east,
|f| f.east(),
eb.west.as_mut(),
time,
),
@ -994,7 +994,7 @@ fn vortexify(
fiter.next().unwrap(),
);
let (y, x) = yx;
vparams.evaluate(time, x, y, rho, rhou, rhov, e);
vparams.evaluate(time, x, y, rho, rhou, rhov, e)
}
#[allow(non_snake_case)]

10
gridgeneration/Cargo.toml
View File

@ -2,13 +2,13 @@
name = "gridgeneration"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
ndarray = { version = "0.15.4", default-features = false }
plotters = { version = "0.3.1", default-features = false, features = ["svg_backend", "line_series", "point_series"] }
ndarray = { version = "0.15.0", default-features = false }
plotters = { version = "0.3.0", default-features = false, features = ["svg_backend", "line_series", "point_series"] }
sbp = { path = "../sbp" }
json5 = { version = "0.4.1", optional = true }
json5 = { version = "0.2.8", optional = true }
[features]
serde = ["sbp/serde1", "dep:json5"]
serde = ["sbp/serde1", "json5"]

8
heat-equation/Cargo.toml
View File

@ -2,14 +2,14 @@
name = "heat-equation"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
sbp = { path = "../sbp", features = ["sparse"] }
ndarray = "0.15.4"
plotters = { version = "0.3.1", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }
sprs = { version = "0.11.0", default-features = false }
ndarray = "0.15.0"
plotters = { version = "0.3.0", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }
sprs = { version = "0.10.0", default-features = false }
integrate = { path = "../utils/integrate" }
[dev-dependencies]

10
maxwell/Cargo.toml
View File

@ -2,19 +2,19 @@
name = "maxwell"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[features]
sparse = ["sbp/sparse", "dep:sprs"]
sparse = ["sbp/sparse", "sprs"]
[dependencies]
ndarray = "0.15.4"
ndarray = "0.15.0"
sbp = { path = "../sbp" }
sprs = { version = "0.11.0", optional = true, default-features = false }
sprs = { version = "0.10.0", optional = true, default-features = false }
integrate = { path = "../utils/integrate" }
[dev-dependencies]
criterion = "0.4.0"
criterion = "0.3.2"
[[bench]]
name = "bench"

24
maxwell/src/lib.rs
View File

@ -15,13 +15,13 @@ impl Clone for Field {
Self(self.0.clone())
}
fn clone_from(&mut self, source: &Self) {
self.0.clone_from(&source.0);
self.0.clone_from(&source.0)
}
}
impl integrate::Integrable for Field {
type State = Self;
type Diff = Self;
type State = Field;
type Diff = Field;
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
s.0.scaled_add(scale, &o.0);
@ -153,12 +153,12 @@ impl<SBP: SbpOperator2d> System<SBP> {
let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
RHS(op, fut, prev, grid, metrics, wb);
};
let mut time = 0.0;
let mut _time = 0.0;
integrate::integrate::<integrate::Rk4, Field, _>(
rhs_adaptor,
&self.sys.0,
&mut self.sys.1,
&mut time,
&mut _time,
dt,
&mut self.wb.k,
);
@ -213,12 +213,12 @@ impl<UO: SbpOperator2d + UpwindOperator2d> System<UO> {
let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
RHS_upwind(op, fut, prev, grid, metrics, wb);
};
let mut time = 0.0;
let mut _time = 0.0;
integrate::integrate::<integrate::Rk4, Field, _>(
rhs_adaptor,
&self.sys.0,
&mut self.sys.1,
&mut time,
&mut _time,
dt,
&mut self.wb.k,
);
@ -487,8 +487,11 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
metrics: &Metrics,
boundaries: &BoundaryTerms,
) {
let ny = y.ny();
let nx = y.nx();
fn positive_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
let r = Float::hypot(kx, ky);
let r = (kx * kx + ky * ky).sqrt();
[
[ky * ky / r / 2.0, ky / 2.0, -kx * ky / r / 2.0],
[ky / 2.0, r / 2.0, -kx / 2.0],
@ -496,7 +499,7 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
]
}
fn negative_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
let r = Float::hypot(kx, ky);
let r = (kx * kx + ky * ky).sqrt();
[
[-ky * ky / r / 2.0, ky / 2.0, kx * ky / r / 2.0],
[ky / 2.0, -r / 2.0, -kx / 2.0],
@ -504,9 +507,6 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
]
}
let ny = y.ny();
let nx = y.nx();
{
let g = match boundaries.east {
Boundary::This => y.slice(s![.., .., 0]),

4
multigrid/.gitignore vendored
View File

@ -1,4 +0,0 @@
build
dist
*.egg-info
.venv

28
multigrid/Cargo.toml
View File

@ -2,7 +2,7 @@
name = "multigrid"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
@ -10,16 +10,16 @@ sbp = { path = "../sbp", features = ["serde1"] }
euler = { path = "../euler", features = ["serde1"] }
hdf5 = "0.8.1"
integrate = { path = "../utils/integrate" }
rayon = "1.5.3"
indicatif = "0.17.0"
ndarray = { version = "0.15.4", features = ["serde"] }
serde = { version = "1.0.138", features = ["derive"] }
json5 = "0.4.1"
indexmap = { version = "1.9.1", features = ["serde-1"] }
argh = "0.1.8"
evalexpr = "7.2.0"
crossbeam-channel = "0.5.5"
crossbeam-utils = "0.8.10"
parking_lot = "0.12.1"
lock_api = "0.4.7"
arrayvec = "0.7.2"
rayon = "1.3.0"
indicatif = "0.17.0-beta.1"
ndarray = { version = "0.15.0", features = ["serde"] }
serde = { version = "1.0.115", features = ["derive"] }
json5 = "0.3.0"
indexmap = { version = "1.5.2", features = ["serde-1"] }
argh = "0.1.4"
evalexpr = "6.3.0"
crossbeam-channel = "0.5.0"
crossbeam-utils = "0.8.5"
parking_lot = "0.11.1"
lock_api = "0.4.4"
arrayvec = "0.7.1"

6
multigrid/plotter/plotter.py → multigrid/eulerplot
View File

@ -267,7 +267,7 @@ def read_from_file(filename):
return grids, file["t"]
def main():
if __name__ == "__main__":
parser = ArgumentParser(description="Plot a solution from the eulersolver")
parser.add_argument("filename", metavar="filename", type=str)
parser.add_argument("-s", help="Save figure", action="store_true", dest="save")
@ -295,7 +295,3 @@ def main():
plot_pressure_slider(grids, t)
else:
plot_pressure(grids, args.save, args.output)
if __name__ == "__main__":
main()

18
multigrid/pyproject.toml
View File

@ -1,18 +0,0 @@
[project]
name = "plotter"
version = "0.0.1"
dependencies = [
"matplotlib",
"numpy",
"h5py",
]
[project.scripts]
eulerplot = "plotter.plotter:main"
[tool.setuptools]
packages = ["plotter"]
[build-system]
requires = ["setuptools"]
build-backend = "setuptools.build_meta"

6
multigrid/src/input.rs
View File

@ -117,10 +117,10 @@ pub enum InterpolationOperator {
NineH2,
}
impl From<InterpolationOperator> for Box<dyn sbp::operators::InterpolationOperator> {
fn from(val: InterpolationOperator) -> Self {
impl Into<Box<dyn sbp::operators::InterpolationOperator>> for InterpolationOperator {
fn into(self) -> Box<dyn sbp::operators::InterpolationOperator> {
use sbp::operators::{Interpolation4, Interpolation8, Interpolation9, Interpolation9h2};
match val {
match self {
InterpolationOperator::Four => Box::new(Interpolation4),
InterpolationOperator::Eight => Box::new(Interpolation8),
InterpolationOperator::Nine => Box::new(Interpolation9),

9
multigrid/src/main.rs
View File

@ -56,13 +56,13 @@ fn main() {
let config: input::Configuration = match json5::from_str(&filecontents) {
Ok(config) => config,
Err(e) => {
eprintln!("Configuration could not be read: {e}");
eprintln!("Configuration could not be read: {}", e);
if let json5::Error::Message {
location: Some(location),
..
} = e
{
eprintln!("\t{location:?}");
eprintln!("\t{:?}", location);
}
return;
}
@ -157,7 +157,7 @@ fn main() {
println!("Time elapsed: {} seconds", duration.as_secs_f64());
}
if let Some(error) = outinfo.error {
println!("Total error: {error:e}");
println!("Total error: {:e}", error);
}
}
}
@ -166,7 +166,6 @@ fn progressbar(ntime: u64) -> indicatif::ProgressBar {
let progressbar = indicatif::ProgressBar::new(ntime);
progressbar.with_style(
indicatif::ProgressStyle::default_bar()
.template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})")
.unwrap(),
.template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})"),
)
}

11
multigrid/src/parsing.rs
View File

@ -139,12 +139,12 @@ impl input::Configuration {
let xi = operators.xi.unwrap_or_else(|| {
default_operators
.xi
.unwrap_or_else(|| panic!("No xi operator found for grid: {name}"))
.unwrap_or_else(|| panic!("No xi operator found for grid: {}", name))
});
let eta = operators.eta.unwrap_or_else(|| {
default_operators
.eta
.unwrap_or_else(|| panic!("No eta operator found for grid: {name}"))
.unwrap_or_else(|| panic!("No eta operator found for grid: {}", name))
});
use input::Operator as op;
@ -183,8 +183,9 @@ impl input::Configuration {
BoundaryConditions::Expressions(expr) => {
euler::BoundaryCharacteristic::Eval(expr.clone() )
}
BoundaryConditions::NotNeeded => panic!(
"Boundary conditions are not available, but needed for grid {name}"
_ => panic!(
"Boundary conditions are not available, but needed for grid {}",
name
),
},
Some(BoundaryType::This) => euler::BoundaryCharacteristic::This,
@ -192,7 +193,7 @@ impl input::Configuration {
if let BoundaryConditions::Vortex(vortex) = &boundary_conditions {
vortex.clone()
} else {
panic!("Wanted vortex boundary conditions not found, needed for grid {name}")
panic!("Wanted vortex boundary conditions not found, needed for grid {}", name)
},
),
Some(BoundaryType::Neighbour(name)) => {

4
multigrid/src/system.rs
View File

@ -200,7 +200,7 @@ impl BaseSystem {
.zip(push)
.enumerate()
{
let builder = std::thread::Builder::new().name(format!("mg: {name}"));
let builder = std::thread::Builder::new().name(format!("mg: {}", name));
let boundary_conditions = bt.map(|bt| match bt {
euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This,
@ -444,7 +444,7 @@ impl System {
for _ in 0..sys.sys.len() {
e += match sys.recv.recv().unwrap() {
(_, MsgToHost::Error(e)) => e,
(_, m) => panic!("Unexpected message: {m:?}"),
(_, m) => panic!("Unexpected message: {:?}", m),
}
}
e

18
sbp/Cargo.toml
View File

@ -2,26 +2,26 @@
name = "sbp"
version = "0.1.1"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
ndarray = "0.15.4"
approx = "0.5.1"
sprs = { version = "0.11.0", optional = true, default-features = false }
serde = { version = "1.0.138", optional = true, default-features = false, features = ["derive"] }
num-traits = "0.2.15"
ndarray = { version = "0.15.0", features = ["approx"] }
approx = "0.4.0"
sprs = { version = "0.10.0", optional = true, default-features = false }
serde = { version = "1.0.115", optional = true, default-features = false, features = ["derive"] }
num-traits = "0.2.14"
float = { path = "../utils/float" }
constmatrix = { path = "../utils/constmatrix" }
core_simd = { git = "https://github.com/rust-lang/portable-simd" }
[features]
# Use f32 as precision, default is f64
f32 = ["float/f32"]
sparse = ["dep:sprs"]
serde1 = ["dep:serde", "ndarray/serde"]
sparse = ["sprs"]
serde1 = ["serde", "ndarray/serde"]
[dev-dependencies]
iai = "0.1.1"
ndarray = { version = "0.15.4", features = ["approx-0_5"] }
[[bench]]
name = "sbpoperators"

16
sbp/src/operators.rs
View File

@ -56,11 +56,11 @@ pub trait SbpOperator2d: Send + Sync {
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
self.op_xi().diff(p, f);
self.op_xi().diff(p, f)
}
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.diffxi(prev.reversed_axes(), fut.reversed_axes());
self.diffxi(prev.reversed_axes(), fut.reversed_axes())
}
fn hxi(&self) -> &'static [Float] {
@ -91,12 +91,12 @@ pub trait UpwindOperator2d: Send + Sync {
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
UpwindOperator2d::op_xi(self).diss(p, f);
UpwindOperator2d::op_xi(self).diss(p, f)
}
}
// Assuming operator is symmetrical x/y
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.dissxi(prev.reversed_axes(), fut.reversed_axes());
self.dissxi(prev.reversed_axes(), fut.reversed_axes())
}
fn op_xi(&self) -> &dyn UpwindOperator1d;
@ -114,10 +114,10 @@ pub trait InterpolationOperator: Send + Sync {
impl SbpOperator2d for (Box<dyn SbpOperator2d>, Box<dyn SbpOperator2d>) {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.1.diffxi(prev, fut);
self.1.diffxi(prev, fut)
}
fn diffeta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.0.diffeta(prev, fut);
self.0.diffeta(prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
@ -136,10 +136,10 @@ impl SbpOperator2d for (Box<dyn SbpOperator2d>, Box<dyn SbpOperator2d>) {
impl UpwindOperator2d for (Box<dyn UpwindOperator2d>, Box<dyn UpwindOperator2d>) {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.1.dissxi(prev, fut);
self.1.dissxi(prev, fut)
}
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
self.0.disseta(prev, fut);
self.0.disseta(prev, fut)
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
self.1.op_xi()

21
sbp/src/operators/algos.rs
View File

@ -2,7 +2,6 @@ use super::*;
use ndarray::s;
use num_traits::Zero;
use std::convert::TryInto;
use std::simd::StdFloat;
pub(crate) use constmatrix::{ColVector, Matrix, RowVector};
@ -112,7 +111,7 @@ pub(crate) fn diff_op_1d_slice<const M: usize, const N: usize, const D: usize>(
a: &Matrix<Float, M, N>,
b: &ColVector<Float, N>,
) {
c.matmul_float_into(a, b);
c.matmul_float_into(a, b)
}
assert_eq!(prev.len(), fut.len());
@ -199,9 +198,9 @@ pub(crate) fn diff_op_1d<const M: usize, const N: usize, const D: usize>(
assert!(nx >= 2 * M);
if let Some((prev, fut)) = prev.as_slice().zip(fut.as_slice_mut()) {
diff_op_1d_slice(matrix, optype, prev, fut);
diff_op_1d_slice(matrix, optype, prev, fut)
} else {
diff_op_1d_fallback(matrix, optype, prev, fut);
diff_op_1d_fallback(matrix, optype, prev, fut)
}
}
@ -259,7 +258,7 @@ pub(crate) fn diff_op_2d_fallback<const M: usize, const N: usize, const D: usize
if d.is_zero() {
continue;
}
fut.scaled_add(idx * d, &id);
fut.scaled_add(idx * d, &id)
}
}
@ -395,9 +394,9 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
let idx = 1.0 / dx;
#[cfg(not(feature = "f32"))]
type SimdT = std::simd::f64x8;
type SimdT = core_simd::f64x8;
#[cfg(feature = "f32")]
type SimdT = std::simd::f32x16;
type SimdT = core_simd::f32x16;
// How many elements that can be simdified
let simdified = SimdT::LANES * (ny / SimdT::LANES);
@ -443,7 +442,7 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
for (iprev, &bl) in bl.iter().enumerate() {
f = index_to_simd(iprev).mul_add(SimdT::splat(bl), f);
}
f *= SimdT::splat(idx);
f *= idx;
fut.clone_from_slice(f.as_array());
}
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
@ -495,7 +494,7 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
let offset = ifut - half_diag_width + id;
f = index_to_simd(offset).mul_add(SimdT::splat(d), f);
}
f *= SimdT::splat(idx);
f *= idx;
fut.clone_from_slice(f.as_array());
}
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
@ -530,7 +529,7 @@ pub(crate) fn diff_op_2d_sliceable<const M: usize, const N: usize, const D: usiz
for (prev, mut fut) in prev.outer_iter().zip(fut.outer_iter_mut()) {
let prev = &prev.as_slice().unwrap()[..nx];
let fut = &mut fut.as_slice_mut().unwrap()[..nx];
diff_op_1d_slice(matrix, optype, prev, fut);
diff_op_1d_slice(matrix, optype, prev, fut)
}
}
@ -548,7 +547,7 @@ pub(crate) fn diff_op_2d<const M: usize, const N: usize, const D: usize>(
([1, _], [1, _])
if prev.as_slice_memory_order().is_some() && fut.as_slice_memory_order().is_some() =>
{
diff_op_2d_sliceable_y_simd(matrix, optype, prev, fut);
diff_op_2d_sliceable_y_simd(matrix, optype, prev, fut)
}
_ => diff_op_2d_fallback(matrix, optype, prev, fut),
}

4
sbp/src/operators/interpolation/interpolation4.rs
View File

@ -79,7 +79,7 @@ impl InterpolationOperator for Interpolation4 {
ndarray::arr2(Self::F2C_BLOCK).view(),
ndarray::arr2(Self::F2C_DIAG).view(),
(3, 2),
);
)
}
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
assert_eq!(fine.len(), 2 * coarse.len() - 1);
@ -89,7 +89,7 @@ impl InterpolationOperator for Interpolation4 {
ndarray::arr2(Self::C2F_BLOCK).view(),
ndarray::arr2(Self::C2F_DIAG).view(),
(4, 1),
);
)
}
}

4
sbp/src/operators/interpolation/interpolation8.rs
View File

@ -58,7 +58,7 @@ impl InterpolationOperator for Interpolation8 {
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
.unwrap();
super::interpolate(fine, coarse, block.view(), diag.view(), (0, 2));
super::interpolate(fine, coarse, block.view(), diag.view(), (0, 2))
}
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
assert_eq!(fine.len(), 2 * coarse.len() - 1);
@ -69,7 +69,7 @@ impl InterpolationOperator for Interpolation8 {
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
.unwrap();
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
}
}

4
sbp/src/operators/interpolation/interpolation9.rs
View File

@ -59,7 +59,7 @@ impl InterpolationOperator for Interpolation9 {
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
.unwrap();
super::interpolate(fine, coarse, block.view(), diag.view(), (5, 2));
super::interpolate(fine, coarse, block.view(), diag.view(), (5, 2))
}
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
assert_eq!(fine.len(), 2 * coarse.len() - 1);
@ -70,7 +70,7 @@ impl InterpolationOperator for Interpolation9 {
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
.unwrap();
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
}
}

4
sbp/src/operators/interpolation/interpolation9h2.rs
View File

@ -59,7 +59,7 @@ impl InterpolationOperator for Interpolation9h2 {
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
.unwrap();
super::interpolate(fine, coarse, block.view(), diag.view(), (4, 2));
super::interpolate(fine, coarse, block.view(), diag.view(), (4, 2))
}
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
assert_eq!(fine.len(), 2 * (coarse.len() - 1));
@ -70,7 +70,7 @@ impl InterpolationOperator for Interpolation9h2 {
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
.unwrap();
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
}
}

6
sbp/src/operators/traditional4.rs
View File

@ -50,7 +50,7 @@ impl SBP4 {
impl SbpOperator1d for SBP4 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
}
fn h(&self) -> &'static [Float] {
@ -74,7 +74,7 @@ impl SbpOperator1d for SBP4 {
impl SbpOperator2d for SBP4 {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut);
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
@ -85,7 +85,7 @@ impl super::SbpOperator1d2 for SBP4 {
fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::D2, OperatorType::Normal, prev, fut.view_mut());
let hi = (prev.len() - 1) as Float;
fut.map_inplace(|x| *x *= hi);
fut.map_inplace(|x| *x *= hi)
}
fn d1(&self) -> &[Float] {
&[-88.0 / 17.0, 144.0 / 17.0, -72.0 / 17.0, 16.0 / 17.0]

4
sbp/src/operators/traditional8.rs
View File

@ -36,7 +36,7 @@ impl SBP8 {
impl SbpOperator1d for SBP8 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
}
fn h(&self) -> &'static [Float] {
@ -56,7 +56,7 @@ impl SbpOperator1d for SBP8 {
impl SbpOperator2d for SBP8 {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut);
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self

6
sbp/src/operators/upwind4.rs
View File

@ -51,7 +51,7 @@ impl Upwind4 {
impl SbpOperator1d for Upwind4 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
}
fn h(&self) -> &'static [Float] {
&Self::H.start
@ -173,7 +173,7 @@ fn upwind4_test() {
impl UpwindOperator1d for Upwind4 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut)
}
fn as_sbp(&self) -> &dyn SbpOperator1d {
@ -190,7 +190,7 @@ impl UpwindOperator2d for Upwind4 {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut);
super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut)
}
fn op_xi(&self) -> &dyn UpwindOperator1d {

8
sbp/src/operators/upwind4h2.rs
View File

@ -53,7 +53,7 @@ impl Upwind4h2 {
impl SbpOperator1d for Upwind4h2 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, OperatorType::H2, prev, fut);
super::diff_op_1d(&Self::DIFF, OperatorType::H2, prev, fut)
}
fn h(&self) -> &'static [Float] {
@ -81,7 +81,7 @@ impl SbpOperator2d for Upwind4h2 {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut);
super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
@ -95,7 +95,7 @@ impl UpwindOperator2d for Upwind4h2 {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut);
super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut)
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self
@ -129,7 +129,7 @@ fn upwind4h2_test() {
impl UpwindOperator1d for Upwind4h2 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut);
super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut)
}
fn as_sbp(&self) -> &dyn SbpOperator1d {

8
sbp/src/operators/upwind9.rs
View File

@ -59,7 +59,7 @@ impl Upwind9 {
impl SbpOperator1d for Upwind9 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DIFF, super::OperatorType::Normal, prev, fut)
}
fn h(&self) -> &'static [Float] {
@ -84,7 +84,7 @@ impl SbpOperator2d for Upwind9 {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut);
super::diff_op_2d(&Self::DIFF, OperatorType::Normal, prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
@ -96,7 +96,7 @@ impl SbpOperator2d for Upwind9 {
impl UpwindOperator1d for Upwind9 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut);
super::diff_op_1d(&Self::DISS, super::OperatorType::Normal, prev, fut)
}
fn as_sbp(&self) -> &dyn SbpOperator1d {
@ -112,7 +112,7 @@ impl UpwindOperator1d for Upwind9 {
impl UpwindOperator2d for Upwind9 {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut);
super::diff_op_2d(&Self::DISS, OperatorType::Normal, prev, fut)
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self

8
sbp/src/operators/upwind9h2.rs
View File

@ -59,7 +59,7 @@ impl Upwind9h2 {
impl SbpOperator1d for Upwind9h2 {
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DIFF, super::OperatorType::H2, prev, fut);
super::diff_op_1d(&Self::DIFF, super::OperatorType::H2, prev, fut)
}
fn h(&self) -> &'static [Float] {
@ -87,7 +87,7 @@ impl SbpOperator2d for Upwind9h2 {
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut);
super::diff_op_2d(&Self::DIFF, OperatorType::H2, prev, fut)
}
fn op_xi(&self) -> &dyn SbpOperator1d {
&Self
@ -124,7 +124,7 @@ fn upwind9h2_test() {
impl UpwindOperator1d for Upwind9h2 {
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut);
super::diff_op_1d(&Self::DISS, super::OperatorType::H2, prev, fut)
}
fn as_sbp(&self) -> &dyn SbpOperator1d {
self
@ -140,7 +140,7 @@ impl UpwindOperator2d for Upwind9h2 {
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
assert_eq!(prev.shape(), fut.shape());
super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut);
super::diff_op_2d(&Self::DISS, OperatorType::H2, prev, fut)
}
fn op_xi(&self) -> &dyn UpwindOperator1d {
&Self

4
sbp/src/utils.rs
View File

@ -4,12 +4,12 @@ use crate::Float;
mod jacobi;
#[cfg(feature = "sparse")]
pub use jacobi::*;
#[cfg(feature = "serde1")]
#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};
#[cfg(feature = "sparse")]
pub use sprs::kronecker_product;
#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Copy, Clone, Debug, Default)]
/// struct to hold output for four directions
pub struct Direction<T> {

6
shallow_water/Cargo.toml
View File

@ -2,10 +2,10 @@
name = "shallow-water"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
ndarray = "0.15.4"
ndarray = "0.15.0"
sbp = { path = "../sbp" }
log = "0.4.17"
log = "0.4.8"
integrate = { path = "../utils/integrate" }

4
shallow_water/src/lib.rs
View File

@ -19,8 +19,8 @@ impl Clone for Field {
}
impl integrate::Integrable for Field {
type State = Self;
type Diff = Self;
type State = Field;
type Diff = Field;
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
s.0.scaled_add(scale, &o.0);

6
utils/constmatrix/Cargo.toml
View File

@ -2,10 +2,10 @@
name = "constmatrix"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
approx = { version = "0.5.1", optional = true }
approx = { version = "0.4.0", optional = true }
cfg-if = "1.0.0"
float = { path = "../float" }
num-traits = "0.2.15"
num-traits = "0.2.14"

2
utils/constmatrix/src/lib.rs
View File

@ -258,7 +258,7 @@ mod approx {
fn default_epsilon() -> Self::Epsilon {
T::default_epsilon()
}
fn abs_diff_eq(&self, other: &Self, _epsilon: Self::Epsilon) -> bool {
fn abs_diff_eq(&self, other: &Self, epsilon: Self::Epsilon) -> bool {
self.iter()
.zip(other.iter())
.all(|(r, l)| r.abs_diff_eq(l, T::default_epsilon()))

6
utils/fast-float/Cargo.toml
View File

@ -2,8 +2,10 @@
name = "fast-float"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
[dependencies]
float = { path = "../float" }
num-traits = "0.2.15"
num-traits = "0.2.14"

2
utils/float/Cargo.toml
View File

@ -2,7 +2,7 @@
name = "float"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[features]
f32 = []

2
utils/integrate/Cargo.toml
View File

@ -2,7 +2,7 @@
name = "integrate"
version = "0.1.0"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[dependencies]
float = { path = "../float/" }

11
webfront/Cargo.toml
View File

@ -2,18 +2,19 @@
name = "sbp-web"
version = "0.1.1"
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
edition = "2021"
edition = "2018"
[lib]
crate-type = ["cdylib"]
[dependencies]
wasm-bindgen = "0.2.81"
console_error_panic_hook = "0.1.7"
wasm-bindgen = "0.2.63"
console_error_panic_hook = "0.1.6"
wee_alloc = "0.4.5"
sbp = { path = "../sbp", features = ["f32"] }
ndarray = "0.15.4"
ndarray = "0.15.0"
euler = { path = "../euler" }
maxwell = { path = "../maxwell" }
shallow-water = { path = "../shallow_water" }
console_log = "0.2.0"
log = "0.4.17"
log = "0.4.8"

4
webfront/src/lib.rs
View File

@ -9,6 +9,10 @@ pub use crate::maxwell::MaxwellUniverse;
mod shallow_water;
pub use crate::shallow_water::ShallowWaterUniverse;
#[cfg(feature = "wee_alloc")]
#[global_allocator]
static ALLOC: wee_alloc::WeeAlloc = wee_alloc::WeeAlloc::INIT;
#[wasm_bindgen]
pub fn set_panic_hook() {
console_error_panic_hook::set_once();