Compare commits
10 Commits
cfeb30fac0
...
11807ffb68
Author | SHA1 | Date |
---|---|---|
Magnus Ulimoen | 11807ffb68 | |
Magnus Ulimoen | 9010a47420 | |
Magnus Ulimoen | 909e15572e | |
Magnus Ulimoen | 4c5c0305e4 | |
Magnus Ulimoen | 5acd46af6d | |
Magnus Ulimoen | 2a1bb3f815 | |
Magnus Ulimoen | bb1909c2a8 | |
Magnus Ulimoen | f40de866ce | |
Magnus Ulimoen | 70cab01334 | |
Magnus Ulimoen | 05873c0e22 |
|
@ -2,7 +2,7 @@
|
||||||
name = "euler"
|
name = "euler"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
# Internal feature flag to gate the expensive tests
|
# Internal feature flag to gate the expensive tests
|
||||||
|
@ -11,15 +11,15 @@ expensive_tests = []
|
||||||
serde1 = ["serde", "arrayvec/serde"]
|
serde1 = ["serde", "arrayvec/serde"]
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
ndarray = "0.15.0"
|
ndarray = "0.15.4"
|
||||||
sbp = { path = "../sbp" }
|
sbp = { path = "../sbp" }
|
||||||
arrayvec = "0.6.0"
|
arrayvec = "0.7.2"
|
||||||
serde = { version = "1.0.115", default-features = false, optional = true, features = ["derive"] }
|
serde = { version = "1.0.138", default-features = false, optional = true, features = ["derive"] }
|
||||||
integrate = { path = "../utils/integrate" }
|
integrate = { path = "../utils/integrate" }
|
||||||
once_cell = "1.7.2"
|
once_cell = "1.13.0"
|
||||||
|
|
||||||
[dev-dependencies]
|
[dev-dependencies]
|
||||||
criterion = "0.3.2"
|
criterion = "0.4.0"
|
||||||
|
|
||||||
[[bench]]
|
[[bench]]
|
||||||
name = "bench"
|
name = "bench"
|
||||||
|
|
|
@ -60,7 +60,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
|
||||||
let metrics = &self.grid.1;
|
let metrics = &self.grid.1;
|
||||||
let rhs_trad = |k: &mut Diff, y: &Field, _time: Float| {
|
let rhs_trad = |k: &mut Diff, y: &Field, _time: Float| {
|
||||||
let boundaries = boundary_extractor(y, grid, &bc);
|
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, _>(
|
integrate::integrate::<integrate::Rk4, Field, _>(
|
||||||
rhs_trad,
|
rhs_trad,
|
||||||
|
@ -98,7 +98,7 @@ impl<SBP: SbpOperator2d> System<SBP> {
|
||||||
|
|
||||||
self.sys
|
self.sys
|
||||||
.0
|
.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 {
|
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 rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
|
||||||
let (grid, metrics) = grid;
|
let (grid, metrics) = grid;
|
||||||
let boundaries = boundary_extractor(y, grid, &bc);
|
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, _>(
|
integrate::integrate::<integrate::Rk4, Field, _>(
|
||||||
rhs_upwind,
|
rhs_upwind,
|
||||||
|
@ -159,7 +159,7 @@ impl<UO: UpwindOperator2d + SbpOperator2d> System<UO> {
|
||||||
let mut rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
|
let mut rhs_upwind = |k: &mut Diff, y: &Field, _time: Float| {
|
||||||
let (grid, metrics) = grid;
|
let (grid, metrics) = grid;
|
||||||
let boundaries = boundary_extractor(y, grid, &bc);
|
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 time = 0.0;
|
||||||
let mut sys2 = self.sys.0.clone();
|
let mut sys2 = self.sys.0.clone();
|
||||||
|
@ -194,7 +194,7 @@ impl Clone for Field {
|
||||||
Self(self.0.clone())
|
Self(self.0.clone())
|
||||||
}
|
}
|
||||||
fn clone_from(&mut self, source: &Self) {
|
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,
|
vortex_param: &VortexParameters,
|
||||||
) {
|
) {
|
||||||
let (rho, rhou, rhov, e) = self.components_mut();
|
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)]
|
#[allow(clippy::erasing_op, clippy::identity_op)]
|
||||||
fn iter(&self) -> impl ExactSizeIterator<Item = FieldValue> + '_ {
|
fn iter(&self) -> impl ExactSizeIterator<Item = FieldValue> + '_ {
|
||||||
|
@ -508,7 +508,7 @@ pub fn RHS_trad(
|
||||||
eflux in &dE.0,
|
eflux in &dE.0,
|
||||||
fflux in &dF.0,
|
fflux in &dF.0,
|
||||||
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
|
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);
|
SAT_characteristics(op, k, y, metrics, boundaries);
|
||||||
|
@ -555,14 +555,14 @@ pub fn RHS_no_SAT(
|
||||||
ad_xi in &ad_xi.0,
|
ad_xi in &ad_xi.0,
|
||||||
ad_eta in &ad_eta.0,
|
ad_eta in &ad_eta.0,
|
||||||
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
|
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 {
|
} else {
|
||||||
azip!((out in &mut k.0,
|
azip!((out in &mut k.0,
|
||||||
eflux in &dE.0,
|
eflux in &dE.0,
|
||||||
fflux in &dF.0,
|
fflux in &dF.0,
|
||||||
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
|
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_xi in &ad_xi.0,
|
||||||
ad_eta in &ad_eta.0,
|
ad_eta in &ad_eta.0,
|
||||||
detj in &metrics.detj().broadcast((4, y.ny(), y.nx())).unwrap()) {
|
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);
|
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();
|
let mut p = wb.rho_mut();
|
||||||
azip!((p in &mut p, &rho in &rho, &rhou in &rhou, &rhov in &rhov, &e in &e) {
|
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();
|
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;
|
let fflux = *ff;
|
||||||
*ef = j_dxi_dx * eflux + j_dxi_dy * fflux;
|
*ef = j_dxi_dx * eflux + j_dxi_dy * fflux;
|
||||||
*ff = j_deta_dx * eflux + j_deta_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 {
|
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||||
match self {
|
match self {
|
||||||
Self::This => write!(f, "This"),
|
Self::This => write!(f, "This"),
|
||||||
Self::Grid(g) => write!(f, "Grid({})", g),
|
Self::Grid(g) => write!(f, "Grid({g})"),
|
||||||
Self::Vortex(vp) => write!(f, "{:?}", vp),
|
Self::Vortex(vp) => write!(f, "{vp:?}"),
|
||||||
Self::Eval(_) => write!(f, "Eval"),
|
Self::Eval(_) => write!(f, "Eval"),
|
||||||
Self::Interpolate(_, _) => write!(f, "Interpolate"),
|
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])
|
to.slice_mut(s![.., i..i + n])
|
||||||
.assign(&seldir(&fields[g]).slice(s![.., start..end]));
|
.assign(&seldir(&fields[g]).slice(s![.., start..end]));
|
||||||
remaining -= 1;
|
remaining -= 1;
|
||||||
if remaining != 0 {
|
if remaining == 0 {
|
||||||
to.slice_mut(s![.., i]).iter_mut().for_each(|x| *x /= 2.0);
|
|
||||||
i += n - 1;
|
|
||||||
} else {
|
|
||||||
i += n;
|
i += n;
|
||||||
assert_eq!(i, to.len_of(Axis(1)));
|
assert_eq!(i, to.len_of(Axis(1)));
|
||||||
|
} else {
|
||||||
|
to.slice_mut(s![.., i]).iter_mut().for_each(|x| *x /= 2.0);
|
||||||
|
i += n - 1;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
to.view()
|
to.view()
|
||||||
|
@ -880,7 +880,7 @@ pub fn boundary_extracts<'a>(
|
||||||
&bt.north,
|
&bt.north,
|
||||||
field,
|
field,
|
||||||
grid.north(),
|
grid.north(),
|
||||||
|f| f.south(),
|
Field::south,
|
||||||
eb.north.as_mut(),
|
eb.north.as_mut(),
|
||||||
time,
|
time,
|
||||||
),
|
),
|
||||||
|
@ -889,7 +889,7 @@ pub fn boundary_extracts<'a>(
|
||||||
&bt.south,
|
&bt.south,
|
||||||
field,
|
field,
|
||||||
grid.south(),
|
grid.south(),
|
||||||
|f| f.north(),
|
Field::north,
|
||||||
eb.south.as_mut(),
|
eb.south.as_mut(),
|
||||||
time,
|
time,
|
||||||
),
|
),
|
||||||
|
@ -898,7 +898,7 @@ pub fn boundary_extracts<'a>(
|
||||||
&bt.east,
|
&bt.east,
|
||||||
field,
|
field,
|
||||||
grid.east(),
|
grid.east(),
|
||||||
|f| f.west(),
|
Field::west,
|
||||||
eb.east.as_mut(),
|
eb.east.as_mut(),
|
||||||
time,
|
time,
|
||||||
),
|
),
|
||||||
|
@ -907,7 +907,7 @@ pub fn boundary_extracts<'a>(
|
||||||
&bt.west,
|
&bt.west,
|
||||||
field,
|
field,
|
||||||
grid.west(),
|
grid.west(),
|
||||||
|f| f.east(),
|
Field::east,
|
||||||
eb.west.as_mut(),
|
eb.west.as_mut(),
|
||||||
time,
|
time,
|
||||||
),
|
),
|
||||||
|
@ -994,7 +994,7 @@ fn vortexify(
|
||||||
fiter.next().unwrap(),
|
fiter.next().unwrap(),
|
||||||
);
|
);
|
||||||
let (y, x) = yx;
|
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)]
|
#[allow(non_snake_case)]
|
||||||
|
|
|
@ -2,13 +2,13 @@
|
||||||
name = "gridgeneration"
|
name = "gridgeneration"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
ndarray = { version = "0.15.0", default-features = false }
|
ndarray = { version = "0.15.4", default-features = false }
|
||||||
plotters = { version = "0.3.0", default-features = false, features = ["svg_backend", "line_series", "point_series"] }
|
plotters = { version = "0.3.1", default-features = false, features = ["svg_backend", "line_series", "point_series"] }
|
||||||
sbp = { path = "../sbp" }
|
sbp = { path = "../sbp" }
|
||||||
json5 = { version = "0.2.8", optional = true }
|
json5 = { version = "0.4.1", optional = true }
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
serde = ["sbp/serde1", "json5"]
|
serde = ["sbp/serde1", "dep:json5"]
|
||||||
|
|
|
@ -2,14 +2,14 @@
|
||||||
name = "heat-equation"
|
name = "heat-equation"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
sbp = { path = "../sbp", features = ["sparse"] }
|
sbp = { path = "../sbp", features = ["sparse"] }
|
||||||
ndarray = "0.15.0"
|
ndarray = "0.15.4"
|
||||||
plotters = { version = "0.3.0", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }
|
plotters = { version = "0.3.1", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] }
|
||||||
sprs = { version = "0.10.0", default-features = false }
|
sprs = { version = "0.11.0", default-features = false }
|
||||||
integrate = { path = "../utils/integrate" }
|
integrate = { path = "../utils/integrate" }
|
||||||
|
|
||||||
[dev-dependencies]
|
[dev-dependencies]
|
||||||
|
|
|
@ -2,19 +2,19 @@
|
||||||
name = "maxwell"
|
name = "maxwell"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
sparse = ["sbp/sparse", "sprs"]
|
sparse = ["sbp/sparse", "dep:sprs"]
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
ndarray = "0.15.0"
|
ndarray = "0.15.4"
|
||||||
sbp = { path = "../sbp" }
|
sbp = { path = "../sbp" }
|
||||||
sprs = { version = "0.10.0", optional = true, default-features = false }
|
sprs = { version = "0.11.0", optional = true, default-features = false }
|
||||||
integrate = { path = "../utils/integrate" }
|
integrate = { path = "../utils/integrate" }
|
||||||
|
|
||||||
[dev-dependencies]
|
[dev-dependencies]
|
||||||
criterion = "0.3.2"
|
criterion = "0.4.0"
|
||||||
|
|
||||||
[[bench]]
|
[[bench]]
|
||||||
name = "bench"
|
name = "bench"
|
||||||
|
|
|
@ -15,13 +15,13 @@ impl Clone for Field {
|
||||||
Self(self.0.clone())
|
Self(self.0.clone())
|
||||||
}
|
}
|
||||||
fn clone_from(&mut self, source: &Self) {
|
fn clone_from(&mut self, source: &Self) {
|
||||||
self.0.clone_from(&source.0)
|
self.0.clone_from(&source.0);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
impl integrate::Integrable for Field {
|
impl integrate::Integrable for Field {
|
||||||
type State = Field;
|
type State = Self;
|
||||||
type Diff = Field;
|
type Diff = Self;
|
||||||
|
|
||||||
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
|
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
|
||||||
s.0.scaled_add(scale, &o.0);
|
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| {
|
let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
|
||||||
RHS(op, fut, prev, grid, metrics, wb);
|
RHS(op, fut, prev, grid, metrics, wb);
|
||||||
};
|
};
|
||||||
let mut _time = 0.0;
|
let mut time = 0.0;
|
||||||
integrate::integrate::<integrate::Rk4, Field, _>(
|
integrate::integrate::<integrate::Rk4, Field, _>(
|
||||||
rhs_adaptor,
|
rhs_adaptor,
|
||||||
&self.sys.0,
|
&self.sys.0,
|
||||||
&mut self.sys.1,
|
&mut self.sys.1,
|
||||||
&mut _time,
|
&mut time,
|
||||||
dt,
|
dt,
|
||||||
&mut self.wb.k,
|
&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| {
|
let rhs_adaptor = move |fut: &mut Field, prev: &Field, _time: Float| {
|
||||||
RHS_upwind(op, fut, prev, grid, metrics, wb);
|
RHS_upwind(op, fut, prev, grid, metrics, wb);
|
||||||
};
|
};
|
||||||
let mut _time = 0.0;
|
let mut time = 0.0;
|
||||||
integrate::integrate::<integrate::Rk4, Field, _>(
|
integrate::integrate::<integrate::Rk4, Field, _>(
|
||||||
rhs_adaptor,
|
rhs_adaptor,
|
||||||
&self.sys.0,
|
&self.sys.0,
|
||||||
&mut self.sys.1,
|
&mut self.sys.1,
|
||||||
&mut _time,
|
&mut time,
|
||||||
dt,
|
dt,
|
||||||
&mut self.wb.k,
|
&mut self.wb.k,
|
||||||
);
|
);
|
||||||
|
@ -487,11 +487,8 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
|
||||||
metrics: &Metrics,
|
metrics: &Metrics,
|
||||||
boundaries: &BoundaryTerms,
|
boundaries: &BoundaryTerms,
|
||||||
) {
|
) {
|
||||||
let ny = y.ny();
|
|
||||||
let nx = y.nx();
|
|
||||||
|
|
||||||
fn positive_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
|
fn positive_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
|
||||||
let r = (kx * kx + ky * ky).sqrt();
|
let r = Float::hypot(kx, ky);
|
||||||
[
|
[
|
||||||
[ky * ky / r / 2.0, ky / 2.0, -kx * ky / r / 2.0],
|
[ky * ky / r / 2.0, ky / 2.0, -kx * ky / r / 2.0],
|
||||||
[ky / 2.0, r / 2.0, -kx / 2.0],
|
[ky / 2.0, r / 2.0, -kx / 2.0],
|
||||||
|
@ -499,7 +496,7 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
|
||||||
]
|
]
|
||||||
}
|
}
|
||||||
fn negative_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
|
fn negative_flux(kx: Float, ky: Float) -> [[Float; 3]; 3] {
|
||||||
let r = (kx * kx + ky * ky).sqrt();
|
let r = Float::hypot(kx, ky);
|
||||||
[
|
[
|
||||||
[-ky * ky / r / 2.0, ky / 2.0, kx * ky / r / 2.0],
|
[-ky * ky / r / 2.0, ky / 2.0, kx * ky / r / 2.0],
|
||||||
[ky / 2.0, -r / 2.0, -kx / 2.0],
|
[ky / 2.0, -r / 2.0, -kx / 2.0],
|
||||||
|
@ -507,6 +504,9 @@ fn SAT_characteristics<SBP: SbpOperator2d>(
|
||||||
]
|
]
|
||||||
}
|
}
|
||||||
|
|
||||||
|
let ny = y.ny();
|
||||||
|
let nx = y.nx();
|
||||||
|
|
||||||
{
|
{
|
||||||
let g = match boundaries.east {
|
let g = match boundaries.east {
|
||||||
Boundary::This => y.slice(s![.., .., 0]),
|
Boundary::This => y.slice(s![.., .., 0]),
|
||||||
|
|
|
@ -0,0 +1,4 @@
|
||||||
|
build
|
||||||
|
dist
|
||||||
|
*.egg-info
|
||||||
|
.venv
|
|
@ -2,7 +2,7 @@
|
||||||
name = "multigrid"
|
name = "multigrid"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
|
@ -10,16 +10,16 @@ sbp = { path = "../sbp", features = ["serde1"] }
|
||||||
euler = { path = "../euler", features = ["serde1"] }
|
euler = { path = "../euler", features = ["serde1"] }
|
||||||
hdf5 = "0.8.1"
|
hdf5 = "0.8.1"
|
||||||
integrate = { path = "../utils/integrate" }
|
integrate = { path = "../utils/integrate" }
|
||||||
rayon = "1.3.0"
|
rayon = "1.5.3"
|
||||||
indicatif = "0.17.0-beta.1"
|
indicatif = "0.17.0"
|
||||||
ndarray = { version = "0.15.0", features = ["serde"] }
|
ndarray = { version = "0.15.4", features = ["serde"] }
|
||||||
serde = { version = "1.0.115", features = ["derive"] }
|
serde = { version = "1.0.138", features = ["derive"] }
|
||||||
json5 = "0.3.0"
|
json5 = "0.4.1"
|
||||||
indexmap = { version = "1.5.2", features = ["serde-1"] }
|
indexmap = { version = "1.9.1", features = ["serde-1"] }
|
||||||
argh = "0.1.4"
|
argh = "0.1.8"
|
||||||
evalexpr = "6.3.0"
|
evalexpr = "7.2.0"
|
||||||
crossbeam-channel = "0.5.0"
|
crossbeam-channel = "0.5.5"
|
||||||
crossbeam-utils = "0.8.5"
|
crossbeam-utils = "0.8.10"
|
||||||
parking_lot = "0.11.1"
|
parking_lot = "0.12.1"
|
||||||
lock_api = "0.4.4"
|
lock_api = "0.4.7"
|
||||||
arrayvec = "0.7.1"
|
arrayvec = "0.7.2"
|
||||||
|
|
|
@ -267,7 +267,7 @@ def read_from_file(filename):
|
||||||
return grids, file["t"]
|
return grids, file["t"]
|
||||||
|
|
||||||
|
|
||||||
if __name__ == "__main__":
|
def main():
|
||||||
parser = ArgumentParser(description="Plot a solution from the eulersolver")
|
parser = ArgumentParser(description="Plot a solution from the eulersolver")
|
||||||
parser.add_argument("filename", metavar="filename", type=str)
|
parser.add_argument("filename", metavar="filename", type=str)
|
||||||
parser.add_argument("-s", help="Save figure", action="store_true", dest="save")
|
parser.add_argument("-s", help="Save figure", action="store_true", dest="save")
|
||||||
|
@ -295,3 +295,7 @@ if __name__ == "__main__":
|
||||||
plot_pressure_slider(grids, t)
|
plot_pressure_slider(grids, t)
|
||||||
else:
|
else:
|
||||||
plot_pressure(grids, args.save, args.output)
|
plot_pressure(grids, args.save, args.output)
|
||||||
|
|
||||||
|
|
||||||
|
if __name__ == "__main__":
|
||||||
|
main()
|
|
@ -0,0 +1,18 @@
|
||||||
|
[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"
|
|
@ -117,10 +117,10 @@ pub enum InterpolationOperator {
|
||||||
NineH2,
|
NineH2,
|
||||||
}
|
}
|
||||||
|
|
||||||
impl Into<Box<dyn sbp::operators::InterpolationOperator>> for InterpolationOperator {
|
impl From<InterpolationOperator> for Box<dyn sbp::operators::InterpolationOperator> {
|
||||||
fn into(self) -> Box<dyn sbp::operators::InterpolationOperator> {
|
fn from(val: InterpolationOperator) -> Self {
|
||||||
use sbp::operators::{Interpolation4, Interpolation8, Interpolation9, Interpolation9h2};
|
use sbp::operators::{Interpolation4, Interpolation8, Interpolation9, Interpolation9h2};
|
||||||
match self {
|
match val {
|
||||||
InterpolationOperator::Four => Box::new(Interpolation4),
|
InterpolationOperator::Four => Box::new(Interpolation4),
|
||||||
InterpolationOperator::Eight => Box::new(Interpolation8),
|
InterpolationOperator::Eight => Box::new(Interpolation8),
|
||||||
InterpolationOperator::Nine => Box::new(Interpolation9),
|
InterpolationOperator::Nine => Box::new(Interpolation9),
|
||||||
|
|
|
@ -56,13 +56,13 @@ fn main() {
|
||||||
let config: input::Configuration = match json5::from_str(&filecontents) {
|
let config: input::Configuration = match json5::from_str(&filecontents) {
|
||||||
Ok(config) => config,
|
Ok(config) => config,
|
||||||
Err(e) => {
|
Err(e) => {
|
||||||
eprintln!("Configuration could not be read: {}", e);
|
eprintln!("Configuration could not be read: {e}");
|
||||||
if let json5::Error::Message {
|
if let json5::Error::Message {
|
||||||
location: Some(location),
|
location: Some(location),
|
||||||
..
|
..
|
||||||
} = e
|
} = e
|
||||||
{
|
{
|
||||||
eprintln!("\t{:?}", location);
|
eprintln!("\t{location:?}");
|
||||||
}
|
}
|
||||||
return;
|
return;
|
||||||
}
|
}
|
||||||
|
@ -157,7 +157,7 @@ fn main() {
|
||||||
println!("Time elapsed: {} seconds", duration.as_secs_f64());
|
println!("Time elapsed: {} seconds", duration.as_secs_f64());
|
||||||
}
|
}
|
||||||
if let Some(error) = outinfo.error {
|
if let Some(error) = outinfo.error {
|
||||||
println!("Total error: {:e}", error);
|
println!("Total error: {error:e}");
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -166,6 +166,7 @@ fn progressbar(ntime: u64) -> indicatif::ProgressBar {
|
||||||
let progressbar = indicatif::ProgressBar::new(ntime);
|
let progressbar = indicatif::ProgressBar::new(ntime);
|
||||||
progressbar.with_style(
|
progressbar.with_style(
|
||||||
indicatif::ProgressStyle::default_bar()
|
indicatif::ProgressStyle::default_bar()
|
||||||
.template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})"),
|
.template("{wide_bar:.cyan/blue} {pos}/{len} ({eta})")
|
||||||
|
.unwrap(),
|
||||||
)
|
)
|
||||||
}
|
}
|
||||||
|
|
|
@ -139,12 +139,12 @@ impl input::Configuration {
|
||||||
let xi = operators.xi.unwrap_or_else(|| {
|
let xi = operators.xi.unwrap_or_else(|| {
|
||||||
default_operators
|
default_operators
|
||||||
.xi
|
.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(|| {
|
let eta = operators.eta.unwrap_or_else(|| {
|
||||||
default_operators
|
default_operators
|
||||||
.eta
|
.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;
|
use input::Operator as op;
|
||||||
|
@ -183,9 +183,8 @@ impl input::Configuration {
|
||||||
BoundaryConditions::Expressions(expr) => {
|
BoundaryConditions::Expressions(expr) => {
|
||||||
euler::BoundaryCharacteristic::Eval(expr.clone() )
|
euler::BoundaryCharacteristic::Eval(expr.clone() )
|
||||||
}
|
}
|
||||||
_ => panic!(
|
BoundaryConditions::NotNeeded => panic!(
|
||||||
"Boundary conditions are not available, but needed for grid {}",
|
"Boundary conditions are not available, but needed for grid {name}"
|
||||||
name
|
|
||||||
),
|
),
|
||||||
},
|
},
|
||||||
Some(BoundaryType::This) => euler::BoundaryCharacteristic::This,
|
Some(BoundaryType::This) => euler::BoundaryCharacteristic::This,
|
||||||
|
@ -193,7 +192,7 @@ impl input::Configuration {
|
||||||
if let BoundaryConditions::Vortex(vortex) = &boundary_conditions {
|
if let BoundaryConditions::Vortex(vortex) = &boundary_conditions {
|
||||||
vortex.clone()
|
vortex.clone()
|
||||||
} else {
|
} 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)) => {
|
Some(BoundaryType::Neighbour(name)) => {
|
||||||
|
|
|
@ -200,7 +200,7 @@ impl BaseSystem {
|
||||||
.zip(push)
|
.zip(push)
|
||||||
.enumerate()
|
.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 {
|
let boundary_conditions = bt.map(|bt| match bt {
|
||||||
euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This,
|
euler::BoundaryCharacteristic::This => DistributedBoundaryConditions::This,
|
||||||
|
@ -444,7 +444,7 @@ impl System {
|
||||||
for _ in 0..sys.sys.len() {
|
for _ in 0..sys.sys.len() {
|
||||||
e += match sys.recv.recv().unwrap() {
|
e += match sys.recv.recv().unwrap() {
|
||||||
(_, MsgToHost::Error(e)) => e,
|
(_, MsgToHost::Error(e)) => e,
|
||||||
(_, m) => panic!("Unexpected message: {:?}", m),
|
(_, m) => panic!("Unexpected message: {m:?}"),
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
e
|
e
|
||||||
|
|
|
@ -2,26 +2,26 @@
|
||||||
name = "sbp"
|
name = "sbp"
|
||||||
version = "0.1.1"
|
version = "0.1.1"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
ndarray = { version = "0.15.0", features = ["approx"] }
|
ndarray = "0.15.4"
|
||||||
approx = "0.4.0"
|
approx = "0.5.1"
|
||||||
sprs = { version = "0.10.0", optional = true, default-features = false }
|
sprs = { version = "0.11.0", optional = true, default-features = false }
|
||||||
serde = { version = "1.0.115", optional = true, default-features = false, features = ["derive"] }
|
serde = { version = "1.0.138", optional = true, default-features = false, features = ["derive"] }
|
||||||
num-traits = "0.2.14"
|
num-traits = "0.2.15"
|
||||||
float = { path = "../utils/float" }
|
float = { path = "../utils/float" }
|
||||||
constmatrix = { path = "../utils/constmatrix" }
|
constmatrix = { path = "../utils/constmatrix" }
|
||||||
core_simd = { git = "https://github.com/rust-lang/portable-simd" }
|
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
# Use f32 as precision, default is f64
|
# Use f32 as precision, default is f64
|
||||||
f32 = ["float/f32"]
|
f32 = ["float/f32"]
|
||||||
sparse = ["sprs"]
|
sparse = ["dep:sprs"]
|
||||||
serde1 = ["serde", "ndarray/serde"]
|
serde1 = ["dep:serde", "ndarray/serde"]
|
||||||
|
|
||||||
[dev-dependencies]
|
[dev-dependencies]
|
||||||
iai = "0.1.1"
|
iai = "0.1.1"
|
||||||
|
ndarray = { version = "0.15.4", features = ["approx-0_5"] }
|
||||||
|
|
||||||
[[bench]]
|
[[bench]]
|
||||||
name = "sbpoperators"
|
name = "sbpoperators"
|
||||||
|
|
|
@ -56,11 +56,11 @@ pub trait SbpOperator2d: Send + Sync {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
assert_eq!(prev.shape(), fut.shape());
|
||||||
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
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>) {
|
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] {
|
fn hxi(&self) -> &'static [Float] {
|
||||||
|
@ -91,12 +91,12 @@ pub trait UpwindOperator2d: Send + Sync {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
|
fn dissxi(&self, prev: ArrayView2<Float>, mut fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
assert_eq!(prev.shape(), fut.shape());
|
||||||
for (p, f) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
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
|
// Assuming operator is symmetrical x/y
|
||||||
fn disseta(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
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;
|
fn op_xi(&self) -> &dyn UpwindOperator1d;
|
||||||
|
@ -114,10 +114,10 @@ pub trait InterpolationOperator: Send + Sync {
|
||||||
|
|
||||||
impl SbpOperator2d for (Box<dyn SbpOperator2d>, Box<dyn SbpOperator2d>) {
|
impl SbpOperator2d for (Box<dyn SbpOperator2d>, Box<dyn SbpOperator2d>) {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
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>) {
|
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 {
|
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>) {
|
impl UpwindOperator2d for (Box<dyn UpwindOperator2d>, Box<dyn UpwindOperator2d>) {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
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>) {
|
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 {
|
fn op_xi(&self) -> &dyn UpwindOperator1d {
|
||||||
self.1.op_xi()
|
self.1.op_xi()
|
||||||
|
|
|
@ -2,6 +2,7 @@ use super::*;
|
||||||
use ndarray::s;
|
use ndarray::s;
|
||||||
use num_traits::Zero;
|
use num_traits::Zero;
|
||||||
use std::convert::TryInto;
|
use std::convert::TryInto;
|
||||||
|
use std::simd::StdFloat;
|
||||||
|
|
||||||
pub(crate) use constmatrix::{ColVector, Matrix, RowVector};
|
pub(crate) use constmatrix::{ColVector, Matrix, RowVector};
|
||||||
|
|
||||||
|
@ -111,7 +112,7 @@ pub(crate) fn diff_op_1d_slice<const M: usize, const N: usize, const D: usize>(
|
||||||
a: &Matrix<Float, M, N>,
|
a: &Matrix<Float, M, N>,
|
||||||
b: &ColVector<Float, N>,
|
b: &ColVector<Float, N>,
|
||||||
) {
|
) {
|
||||||
c.matmul_float_into(a, b)
|
c.matmul_float_into(a, b);
|
||||||
}
|
}
|
||||||
|
|
||||||
assert_eq!(prev.len(), fut.len());
|
assert_eq!(prev.len(), fut.len());
|
||||||
|
@ -198,9 +199,9 @@ pub(crate) fn diff_op_1d<const M: usize, const N: usize, const D: usize>(
|
||||||
assert!(nx >= 2 * M);
|
assert!(nx >= 2 * M);
|
||||||
|
|
||||||
if let Some((prev, fut)) = prev.as_slice().zip(fut.as_slice_mut()) {
|
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 {
|
} else {
|
||||||
diff_op_1d_fallback(matrix, optype, prev, fut)
|
diff_op_1d_fallback(matrix, optype, prev, fut);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -258,7 +259,7 @@ pub(crate) fn diff_op_2d_fallback<const M: usize, const N: usize, const D: usize
|
||||||
if d.is_zero() {
|
if d.is_zero() {
|
||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
fut.scaled_add(idx * d, &id)
|
fut.scaled_add(idx * d, &id);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -394,9 +395,9 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
|
||||||
let idx = 1.0 / dx;
|
let idx = 1.0 / dx;
|
||||||
|
|
||||||
#[cfg(not(feature = "f32"))]
|
#[cfg(not(feature = "f32"))]
|
||||||
type SimdT = core_simd::f64x8;
|
type SimdT = std::simd::f64x8;
|
||||||
#[cfg(feature = "f32")]
|
#[cfg(feature = "f32")]
|
||||||
type SimdT = core_simd::f32x16;
|
type SimdT = std::simd::f32x16;
|
||||||
|
|
||||||
// How many elements that can be simdified
|
// How many elements that can be simdified
|
||||||
let simdified = SimdT::LANES * (ny / SimdT::LANES);
|
let simdified = SimdT::LANES * (ny / SimdT::LANES);
|
||||||
|
@ -442,7 +443,7 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
|
||||||
for (iprev, &bl) in bl.iter().enumerate() {
|
for (iprev, &bl) in bl.iter().enumerate() {
|
||||||
f = index_to_simd(iprev).mul_add(SimdT::splat(bl), f);
|
f = index_to_simd(iprev).mul_add(SimdT::splat(bl), f);
|
||||||
}
|
}
|
||||||
f *= idx;
|
f *= SimdT::splat(idx);
|
||||||
fut.clone_from_slice(f.as_array());
|
fut.clone_from_slice(f.as_array());
|
||||||
}
|
}
|
||||||
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
|
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
|
||||||
|
@ -494,7 +495,7 @@ pub(crate) fn diff_op_2d_sliceable_y_simd<const M: usize, const N: usize, const
|
||||||
let offset = ifut - half_diag_width + id;
|
let offset = ifut - half_diag_width + id;
|
||||||
f = index_to_simd(offset).mul_add(SimdT::splat(d), f);
|
f = index_to_simd(offset).mul_add(SimdT::splat(d), f);
|
||||||
}
|
}
|
||||||
f *= idx;
|
f *= SimdT::splat(idx);
|
||||||
fut.clone_from_slice(f.as_array());
|
fut.clone_from_slice(f.as_array());
|
||||||
}
|
}
|
||||||
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
|
for (j, fut) in (simdified..ny).zip(fut.into_remainder()) {
|
||||||
|
@ -529,7 +530,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()) {
|
for (prev, mut fut) in prev.outer_iter().zip(fut.outer_iter_mut()) {
|
||||||
let prev = &prev.as_slice().unwrap()[..nx];
|
let prev = &prev.as_slice().unwrap()[..nx];
|
||||||
let fut = &mut fut.as_slice_mut().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);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -547,7 +548,7 @@ pub(crate) fn diff_op_2d<const M: usize, const N: usize, const D: usize>(
|
||||||
([1, _], [1, _])
|
([1, _], [1, _])
|
||||||
if prev.as_slice_memory_order().is_some() && fut.as_slice_memory_order().is_some() =>
|
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),
|
_ => diff_op_2d_fallback(matrix, optype, prev, fut),
|
||||||
}
|
}
|
||||||
|
|
|
@ -79,7 +79,7 @@ impl InterpolationOperator for Interpolation4 {
|
||||||
ndarray::arr2(Self::F2C_BLOCK).view(),
|
ndarray::arr2(Self::F2C_BLOCK).view(),
|
||||||
ndarray::arr2(Self::F2C_DIAG).view(),
|
ndarray::arr2(Self::F2C_DIAG).view(),
|
||||||
(3, 2),
|
(3, 2),
|
||||||
)
|
);
|
||||||
}
|
}
|
||||||
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
|
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
|
||||||
assert_eq!(fine.len(), 2 * coarse.len() - 1);
|
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_BLOCK).view(),
|
||||||
ndarray::arr2(Self::C2F_DIAG).view(),
|
ndarray::arr2(Self::C2F_DIAG).view(),
|
||||||
(4, 1),
|
(4, 1),
|
||||||
)
|
);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -58,7 +58,7 @@ impl InterpolationOperator for Interpolation8 {
|
||||||
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
||||||
.unwrap();
|
.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>) {
|
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
|
||||||
assert_eq!(fine.len(), 2 * coarse.len() - 1);
|
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())
|
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
||||||
.unwrap();
|
.unwrap();
|
||||||
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
|
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -59,7 +59,7 @@ impl InterpolationOperator for Interpolation9 {
|
||||||
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
||||||
.unwrap();
|
.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>) {
|
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
|
||||||
assert_eq!(fine.len(), 2 * coarse.len() - 1);
|
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())
|
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
||||||
.unwrap();
|
.unwrap();
|
||||||
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
|
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -59,7 +59,7 @@ impl InterpolationOperator for Interpolation9h2 {
|
||||||
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
let diag = Array::from_iter(Self::F2C_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
.into_shape((Self::F2C_DIAG.len(), Self::F2C_DIAG[0].len()))
|
||||||
.unwrap();
|
.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>) {
|
fn coarse2fine(&self, coarse: ArrayView1<Float>, fine: ArrayViewMut1<Float>) {
|
||||||
assert_eq!(fine.len(), 2 * (coarse.len() - 1));
|
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())
|
let diag = Array::from_iter(Self::C2F_DIAG.iter().flatten().copied())
|
||||||
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
.into_shape((Self::C2F_DIAG.len(), Self::C2F_DIAG[0].len()))
|
||||||
.unwrap();
|
.unwrap();
|
||||||
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1))
|
super::interpolate(coarse, fine, block.view(), diag.view(), (8, 1));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -50,7 +50,7 @@ impl SBP4 {
|
||||||
|
|
||||||
impl SbpOperator1d for SBP4 {
|
impl SbpOperator1d for SBP4 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
|
@ -74,7 +74,7 @@ impl SbpOperator1d for SBP4 {
|
||||||
impl SbpOperator2d for SBP4 {
|
impl SbpOperator2d for SBP4 {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn SbpOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
@ -85,7 +85,7 @@ impl super::SbpOperator1d2 for SBP4 {
|
||||||
fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) {
|
fn diff2(&self, prev: ArrayView1<Float>, mut fut: ArrayViewMut1<Float>) {
|
||||||
super::diff_op_1d(&Self::D2, OperatorType::Normal, prev, fut.view_mut());
|
super::diff_op_1d(&Self::D2, OperatorType::Normal, prev, fut.view_mut());
|
||||||
let hi = (prev.len() - 1) as Float;
|
let hi = (prev.len() - 1) as Float;
|
||||||
fut.map_inplace(|x| *x *= hi)
|
fut.map_inplace(|x| *x *= hi);
|
||||||
}
|
}
|
||||||
fn d1(&self) -> &[Float] {
|
fn d1(&self) -> &[Float] {
|
||||||
&[-88.0 / 17.0, 144.0 / 17.0, -72.0 / 17.0, 16.0 / 17.0]
|
&[-88.0 / 17.0, 144.0 / 17.0, -72.0 / 17.0, 16.0 / 17.0]
|
||||||
|
|
|
@ -36,7 +36,7 @@ impl SBP8 {
|
||||||
|
|
||||||
impl SbpOperator1d for SBP8 {
|
impl SbpOperator1d for SBP8 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
|
@ -56,7 +56,7 @@ impl SbpOperator1d for SBP8 {
|
||||||
impl SbpOperator2d for SBP8 {
|
impl SbpOperator2d for SBP8 {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn SbpOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
|
|
@ -51,7 +51,7 @@ impl Upwind4 {
|
||||||
|
|
||||||
impl SbpOperator1d for Upwind4 {
|
impl SbpOperator1d for Upwind4 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
&Self::H.start
|
&Self::H.start
|
||||||
|
@ -173,7 +173,7 @@ fn upwind4_test() {
|
||||||
|
|
||||||
impl UpwindOperator1d for Upwind4 {
|
impl UpwindOperator1d for Upwind4 {
|
||||||
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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 {
|
fn as_sbp(&self) -> &dyn SbpOperator1d {
|
||||||
|
@ -190,7 +190,7 @@ impl UpwindOperator2d for Upwind4 {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn UpwindOperator1d {
|
||||||
|
|
|
@ -53,7 +53,7 @@ impl Upwind4h2 {
|
||||||
|
|
||||||
impl SbpOperator1d for Upwind4h2 {
|
impl SbpOperator1d for Upwind4h2 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
|
@ -81,7 +81,7 @@ impl SbpOperator2d for Upwind4h2 {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn SbpOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
@ -95,7 +95,7 @@ impl UpwindOperator2d for Upwind4h2 {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn UpwindOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
@ -129,7 +129,7 @@ fn upwind4h2_test() {
|
||||||
|
|
||||||
impl UpwindOperator1d for Upwind4h2 {
|
impl UpwindOperator1d for Upwind4h2 {
|
||||||
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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 {
|
fn as_sbp(&self) -> &dyn SbpOperator1d {
|
||||||
|
|
|
@ -59,7 +59,7 @@ impl Upwind9 {
|
||||||
|
|
||||||
impl SbpOperator1d for Upwind9 {
|
impl SbpOperator1d for Upwind9 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
|
@ -84,7 +84,7 @@ impl SbpOperator2d for Upwind9 {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn SbpOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
@ -96,7 +96,7 @@ impl SbpOperator2d for Upwind9 {
|
||||||
|
|
||||||
impl UpwindOperator1d for Upwind9 {
|
impl UpwindOperator1d for Upwind9 {
|
||||||
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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 {
|
fn as_sbp(&self) -> &dyn SbpOperator1d {
|
||||||
|
@ -112,7 +112,7 @@ impl UpwindOperator1d for Upwind9 {
|
||||||
impl UpwindOperator2d for Upwind9 {
|
impl UpwindOperator2d for Upwind9 {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn UpwindOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
|
|
@ -59,7 +59,7 @@ impl Upwind9h2 {
|
||||||
|
|
||||||
impl SbpOperator1d for Upwind9h2 {
|
impl SbpOperator1d for Upwind9h2 {
|
||||||
fn diff(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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] {
|
fn h(&self) -> &'static [Float] {
|
||||||
|
@ -87,7 +87,7 @@ impl SbpOperator2d for Upwind9h2 {
|
||||||
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn diffxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn SbpOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
@ -124,7 +124,7 @@ fn upwind9h2_test() {
|
||||||
|
|
||||||
impl UpwindOperator1d for Upwind9h2 {
|
impl UpwindOperator1d for Upwind9h2 {
|
||||||
fn diss(&self, prev: ArrayView1<Float>, fut: ArrayViewMut1<Float>) {
|
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 {
|
fn as_sbp(&self) -> &dyn SbpOperator1d {
|
||||||
self
|
self
|
||||||
|
@ -140,7 +140,7 @@ impl UpwindOperator2d for Upwind9h2 {
|
||||||
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
fn dissxi(&self, prev: ArrayView2<Float>, fut: ArrayViewMut2<Float>) {
|
||||||
assert_eq!(prev.shape(), fut.shape());
|
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 {
|
fn op_xi(&self) -> &dyn UpwindOperator1d {
|
||||||
&Self
|
&Self
|
||||||
|
|
|
@ -4,12 +4,12 @@ use crate::Float;
|
||||||
mod jacobi;
|
mod jacobi;
|
||||||
#[cfg(feature = "sparse")]
|
#[cfg(feature = "sparse")]
|
||||||
pub use jacobi::*;
|
pub use jacobi::*;
|
||||||
#[cfg(feature = "serde")]
|
#[cfg(feature = "serde1")]
|
||||||
use serde::{Deserialize, Serialize};
|
use serde::{Deserialize, Serialize};
|
||||||
#[cfg(feature = "sparse")]
|
#[cfg(feature = "sparse")]
|
||||||
pub use sprs::kronecker_product;
|
pub use sprs::kronecker_product;
|
||||||
|
|
||||||
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
|
#[cfg_attr(feature = "serde1", derive(Serialize, Deserialize))]
|
||||||
#[derive(Copy, Clone, Debug, Default)]
|
#[derive(Copy, Clone, Debug, Default)]
|
||||||
/// struct to hold output for four directions
|
/// struct to hold output for four directions
|
||||||
pub struct Direction<T> {
|
pub struct Direction<T> {
|
||||||
|
|
|
@ -2,10 +2,10 @@
|
||||||
name = "shallow-water"
|
name = "shallow-water"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
ndarray = "0.15.0"
|
ndarray = "0.15.4"
|
||||||
sbp = { path = "../sbp" }
|
sbp = { path = "../sbp" }
|
||||||
log = "0.4.8"
|
log = "0.4.17"
|
||||||
integrate = { path = "../utils/integrate" }
|
integrate = { path = "../utils/integrate" }
|
||||||
|
|
|
@ -19,8 +19,8 @@ impl Clone for Field {
|
||||||
}
|
}
|
||||||
|
|
||||||
impl integrate::Integrable for Field {
|
impl integrate::Integrable for Field {
|
||||||
type State = Field;
|
type State = Self;
|
||||||
type Diff = Field;
|
type Diff = Self;
|
||||||
|
|
||||||
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
|
fn scaled_add(s: &mut Self::State, o: &Self::Diff, scale: Float) {
|
||||||
s.0.scaled_add(scale, &o.0);
|
s.0.scaled_add(scale, &o.0);
|
||||||
|
|
|
@ -2,10 +2,10 @@
|
||||||
name = "constmatrix"
|
name = "constmatrix"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
approx = { version = "0.4.0", optional = true }
|
approx = { version = "0.5.1", optional = true }
|
||||||
cfg-if = "1.0.0"
|
cfg-if = "1.0.0"
|
||||||
float = { path = "../float" }
|
float = { path = "../float" }
|
||||||
num-traits = "0.2.14"
|
num-traits = "0.2.15"
|
||||||
|
|
|
@ -258,7 +258,7 @@ mod approx {
|
||||||
fn default_epsilon() -> Self::Epsilon {
|
fn default_epsilon() -> Self::Epsilon {
|
||||||
T::default_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()
|
self.iter()
|
||||||
.zip(other.iter())
|
.zip(other.iter())
|
||||||
.all(|(r, l)| r.abs_diff_eq(l, T::default_epsilon()))
|
.all(|(r, l)| r.abs_diff_eq(l, T::default_epsilon()))
|
||||||
|
|
|
@ -2,10 +2,8 @@
|
||||||
name = "fast-float"
|
name = "fast-float"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
|
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
float = { path = "../float" }
|
float = { path = "../float" }
|
||||||
num-traits = "0.2.14"
|
num-traits = "0.2.15"
|
||||||
|
|
|
@ -2,7 +2,7 @@
|
||||||
name = "float"
|
name = "float"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
f32 = []
|
f32 = []
|
||||||
|
|
|
@ -2,7 +2,7 @@
|
||||||
name = "integrate"
|
name = "integrate"
|
||||||
version = "0.1.0"
|
version = "0.1.0"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
float = { path = "../float/" }
|
float = { path = "../float/" }
|
||||||
|
|
|
@ -2,19 +2,18 @@
|
||||||
name = "sbp-web"
|
name = "sbp-web"
|
||||||
version = "0.1.1"
|
version = "0.1.1"
|
||||||
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
authors = ["Magnus Ulimoen <magnus@ulimoen.dev>"]
|
||||||
edition = "2018"
|
edition = "2021"
|
||||||
|
|
||||||
[lib]
|
[lib]
|
||||||
crate-type = ["cdylib"]
|
crate-type = ["cdylib"]
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
wasm-bindgen = "0.2.63"
|
wasm-bindgen = "0.2.81"
|
||||||
console_error_panic_hook = "0.1.6"
|
console_error_panic_hook = "0.1.7"
|
||||||
wee_alloc = "0.4.5"
|
|
||||||
sbp = { path = "../sbp", features = ["f32"] }
|
sbp = { path = "../sbp", features = ["f32"] }
|
||||||
ndarray = "0.15.0"
|
ndarray = "0.15.4"
|
||||||
euler = { path = "../euler" }
|
euler = { path = "../euler" }
|
||||||
maxwell = { path = "../maxwell" }
|
maxwell = { path = "../maxwell" }
|
||||||
shallow-water = { path = "../shallow_water" }
|
shallow-water = { path = "../shallow_water" }
|
||||||
console_log = "0.2.0"
|
console_log = "0.2.0"
|
||||||
log = "0.4.8"
|
log = "0.4.17"
|
||||||
|
|
|
@ -9,10 +9,6 @@ pub use crate::maxwell::MaxwellUniverse;
|
||||||
mod shallow_water;
|
mod shallow_water;
|
||||||
pub use crate::shallow_water::ShallowWaterUniverse;
|
pub use crate::shallow_water::ShallowWaterUniverse;
|
||||||
|
|
||||||
#[cfg(feature = "wee_alloc")]
|
|
||||||
#[global_allocator]
|
|
||||||
static ALLOC: wee_alloc::WeeAlloc = wee_alloc::WeeAlloc::INIT;
|
|
||||||
|
|
||||||
#[wasm_bindgen]
|
#[wasm_bindgen]
|
||||||
pub fn set_panic_hook() {
|
pub fn set_panic_hook() {
|
||||||
console_error_panic_hook::set_once();
|
console_error_panic_hook::set_once();
|
||||||
|
|
Loading…
Reference in New Issue