compute dt from CFL condition

This commit is contained in:
Magnus Ulimoen 2020-05-03 20:45:27 +02:00
parent a1acf00c5a
commit 2f98fae1c6
2 changed files with 48 additions and 14 deletions

View File

@ -8,7 +8,7 @@ edition = "2018"
[dependencies]
sbp = { path = "../sbp", features = ["rayon"] }
euler = { path = "../euler" }
hdf5 = { version = "0.6.0", features = ["static"] }
hdf5 = { version = "0.6.0", features = ["static", "gzip"] }
rayon = "1.3.0"
indicatif = "0.14.0"
structopt = "0.3.13"

View File

@ -129,6 +129,52 @@ impl System {
std::mem::swap(&mut self.fnow, &mut self.fnext);
}
/// Suggested maximum dt for this problem
fn max_dt(&self) -> Float {
let c_max = 1.0;
let mut max_dt: Float = Float::INFINITY;
for (field, metrics) in self.fnow.iter().zip(self.metrics.iter()) {
let nx = field.nx();
let ny = field.ny();
let rho = field.rho();
let rhou = field.rhou();
let rhov = field.rhov();
let mut max_u: Float = 0.0;
let mut max_v: Float = 0.0;
for ((((((rho, rhou), rhov), detj_dxi_dx), detj_dxi_dy), detj_deta_dx), detj_deta_dy) in
rho.iter()
.zip(rhou.iter())
.zip(rhov.iter())
.zip(metrics.detj_dxi_dx())
.zip(metrics.detj_dxi_dy())
.zip(metrics.detj_deta_dx())
.zip(metrics.detj_deta_dy())
{
let u = rhou / rho;
let v = rhov / rho;
let uhat: Float = detj_dxi_dx * u + detj_dxi_dy * v;
let vhat: Float = detj_deta_dx * u + detj_deta_dy * v;
max_u = max_u.max(uhat.abs());
max_v = max_v.max(vhat.abs());
}
let dx = 1.0 / nx as Float;
let dy = 1.0 / ny as Float;
let c_dt = Float::max(max_u / dx, max_v / dy);
max_dt = Float::min(max_dt, c_max / c_dt);
}
max_dt
}
}
#[derive(Debug, StructOpt)]
@ -168,19 +214,7 @@ fn main() {
let mut sys = System::new(grids, bt, operators);
sys.vortex(0.0, &vortexparams);
let max_n = {
let max_nx = sys.grids.iter().map(|g| g.nx()).max().unwrap();
let max_ny = sys.grids.iter().map(|g| g.ny()).max().unwrap();
std::cmp::max(max_nx, max_ny)
};
// Add a robust method for determining CFL, use for example the maximum speed of the initial
// field along with \delta x / \delta y
// U_max = max(rhou/u, rhov/v)
// This requires scaling with the determinant to obtain the "true" speed in computational
// space
// CFL = 0.2
// \delta t = CFL min(\delta x, \delta y) / U_max
let dt = 0.02 / (max_n as Float);
let dt = sys.max_dt();
let ntime = (integration_time / dt).round() as u64;