diff --git a/Cargo.lock b/Cargo.lock index 00efad4..82ad3c2 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -41,7 +41,6 @@ dependencies = [ "ndarray", "ndarray-stats", "pdifflib", - "pdifflib_derive", "queues", "serde", "toml", @@ -403,14 +402,6 @@ dependencies = [ "toml", ] -[[package]] -name = "pdifflib_derive" -version = "0.1.0" -dependencies = [ - "quote", - "syn 2.0.118", -] - [[package]] name = "pkg-config" version = "0.3.33" diff --git a/Cargo.toml b/Cargo.toml index f293350..14c4bd0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -5,7 +5,6 @@ edition = "2021" [dependencies] pdifflib={version= '0.1.0',path="../pdifflib"} -pdifflib_derive={version= '0.1.0',path="../pdifflib_derive"} ndarray = "0.15.1" ndarray-stats = "0.5.0" queues = "1.0.2" diff --git a/src/equations.rs b/src/equations.rs index 89bb668..b0b3457 100644 --- a/src/equations.rs +++ b/src/equations.rs @@ -1,6 +1,6 @@ use crate::{H, NX, NY, PEREODIC}; use ndarray::Array2; -use pdifflib::field::{Field2D, Grid}; +use pdifflib::field::Grid; use pdifflib::finit_diff::poisson_relax; use pdifflib::{dx, dx_b, dx_f, dy, dy_b, dy_f, jacobian, lap, mx_b, my_b}; use serde::{Deserialize, Serialize}; @@ -10,7 +10,7 @@ pub struct Params { pub rel: f64, pub rel_c: f64, pub le: f64, - pub pr: f64, + pub sc: f64, pub pe: f64, pub ma: f64, pub time: f64, @@ -19,20 +19,23 @@ pub struct Params { pub struct TemperatureEquation { grid: Grid, delta: Array2, + le: f64, } impl TemperatureEquation { - pub fn new(nx: usize, ny: usize) -> Self { + pub fn new(nx: usize, ny: usize, le: f64) -> Self { Self { grid: Grid::new(nx, ny), delta: Array2::::zeros((nx, ny)), + le, } } pub fn step(&mut self, temp: &mut Array2, psi: &Array2, dt: f64) { let inv_h2 = 1.0 / (H * H); for (i, j) in self.grid.inner_nodes() { - self.delta[[i, j]] = dt * inv_h2 * (lap!(temp, i, j) + jacobian!(psi, temp, i, j)); + self.delta[[i, j]] = + dt * inv_h2 * (self.le * lap!(temp, i, j) + jacobian!(psi, temp, i, j)); } *temp += &self.delta; } @@ -42,16 +45,20 @@ pub struct PhiEquation { grid: Grid, delta: Array2, rel: f64, - pr: f64, + rel_c: f64, + sc: f64, + le: f64, } impl PhiEquation { - pub fn new(nx: usize, ny: usize, rel: f64, pr: f64) -> Self { + pub fn new(nx: usize, ny: usize, rel: f64, rel_c: f64, sc: f64, le: f64) -> Self { Self { grid: Grid::new(nx, ny), delta: Array2::::zeros((nx, ny)), rel, - pr, + rel_c, + sc, + le, } } @@ -60,14 +67,14 @@ impl PhiEquation { phi: &mut Array2, psi: &Array2, temp: &Array2, + conc: &Array2, dt: f64, ) { let inv_h = 1.0 / H; let inv_h2 = inv_h * inv_h; - let pr = self.pr; for (i, j) in self.grid.inner_nodes() { - let archim = self.rel * dx!(temp, i, j); - let tmp = pr * (lap!(phi, i, j) * inv_h2 + archim * inv_h) + let archim = self.le * self.rel * dx!(temp, i, j) + self.rel_c * dx!(conc, i, j); + let tmp = self.sc * (lap!(phi, i, j) * inv_h2 + archim * inv_h) + jacobian!(psi, phi, i, j) * inv_h2; self.delta[[i, j]] = tmp * dt; } @@ -108,77 +115,83 @@ impl ConcentrationEquation { } } - pub fn step(&mut self, conc: &mut Field2D, psi: &Field2D, _temp: &Field2D, dt: f64) { + pub fn step( + &mut self, + conc: &mut Array2, + psi: &Array2, + _temp: &Array2, + dt: f64, + ) { self.calc_v(psi); let grid = self.grid; for (i, k) in grid.x_faces() { - let mut tmp = self.vx[[i, k]] * mx_b!(conc.f, i, k); - tmp += -dx_b!(conc.f, i, k); + let mut tmp = self.vx[[i, k]] * mx_b!(conc, i, k); + tmp += -dx_b!(conc, i, k); self.qew[[i, k]] = tmp / H; } for (i, k) in grid.y_faces() { - let mut tmp = ((-self.vy[[i, k]]) + self.pe * H) * my_b!(conc.f, i, k); - tmp += -dy_b!(conc.f, i, k); + let mut tmp = ((-self.vy[[i, k]]) + self.pe * H) * my_b!(conc, i, k); + tmp += -dy_b!(conc, i, k); self.qsn[[i, k]] = tmp / H; } for (i, k) in grid.inner_nodes() { let q = -(dx_f!(self.qew, i, k) + dy_f!(self.qsn, i, k)); - conc.f[[i, k]] += dt / H * q; + conc[[i, k]] += dt / H * q; } for i in 1..NX - 1 { let k = 0; - conc.f[[i, k]] += dt / H * (-dx_f!(self.qew, i, k) - 2.0 * self.qsn[[i, k + 1]]); + conc[[i, k]] += dt / H * (-dx_f!(self.qew, i, k) - 2.0 * self.qsn[[i, k + 1]]); let k = NY - 1; - conc.f[[i, k]] += dt / H * (-dx_f!(self.qew, i, k) + 2.0 * self.qsn[[i, k]]); + conc[[i, k]] += dt / H * (-dx_f!(self.qew, i, k) + 2.0 * self.qsn[[i, k]]); } if PEREODIC { for k in 0..NY { - conc.f[[0, k]] = conc.f[[NX - 2, k]]; - conc.f[[NX - 1, k]] = conc.f[[1, k]]; + conc[[0, k]] = conc[[NX - 2, k]]; + conc[[NX - 1, k]] = conc[[1, k]]; } } else { for k in 1..NY - 1 { let i = 0; - conc.f[[i, k]] += dt / H * (-dy_f!(self.qsn, i, k) - 2.0 * self.qew[[i + 1, k]]); + conc[[i, k]] += dt / H * (-dy_f!(self.qsn, i, k) - 2.0 * self.qew[[i + 1, k]]); let i = NX - 1; - conc.f[[i, k]] += dt / H * (-dy_f!(self.qsn, i, k) + 2.0 * self.qew[[i, k]]); + conc[[i, k]] += dt / H * (-dy_f!(self.qsn, i, k) + 2.0 * self.qew[[i, k]]); } - conc.f[[0, 0]] += dt / H * 2.0 * (-self.qew[[1, 0]] - self.qsn[[0, 1]]); - conc.f[[0, NY - 1]] += dt / H * 2.0 * (-self.qew[[1, NY - 1]] + self.qsn[[0, NY - 1]]); - conc.f[[NX - 1, 0]] += dt / H * 2.0 * (self.qew[[NX - 1, 0]] - self.qsn[[NX - 1, 1]]); - conc.f[[NX - 1, NY - 1]] += + conc[[0, 0]] += dt / H * 2.0 * (-self.qew[[1, 0]] - self.qsn[[0, 1]]); + conc[[0, NY - 1]] += dt / H * 2.0 * (-self.qew[[1, NY - 1]] + self.qsn[[0, NY - 1]]); + conc[[NX - 1, 0]] += dt / H * 2.0 * (self.qew[[NX - 1, 0]] - self.qsn[[NX - 1, 1]]); + conc[[NX - 1, NY - 1]] += dt / H * 2.0 * (self.qew[[NX - 1, NY - 1]] + self.qsn[[NX - 1, NY - 1]]); } } - fn calc_v(&mut self, psi: &Field2D) { + fn calc_v(&mut self, psi: &Array2) { for i in 1..NX { let mut k = 0; - self.vx[[i, k]] = 0.5 * (psi.f[[i, k + 1]] + psi.f[[i - 1, k + 1]]); + self.vx[[i, k]] = 0.5 * (psi[[i, k + 1]] + psi[[i - 1, k + 1]]); for k in 1..NY - 1 { - self.vx[[i, k]] = 0.5 * (dy!(psi.f, i, k) + dy!(psi.f, i - 1, k)); + self.vx[[i, k]] = 0.5 * (dy!(psi, i, k) + dy!(psi, i - 1, k)); } k = NY - 1; - self.vx[[i, k]] = -0.5 * (psi.f[[i, k - 1]] + psi.f[[i - 1, k - 1]]); + self.vx[[i, k]] = -0.5 * (psi[[i, k - 1]] + psi[[i - 1, k - 1]]); } for i in 1..NX - 1 { for k in 1..NY { - self.vy[[i, k]] = 0.5 * (dx!(psi.f, i, k) + dx!(psi.f, i, k - 1)); + self.vy[[i, k]] = 0.5 * (dx!(psi, i, k) + dx!(psi, i, k - 1)); } } for k in 1..NY { let i = 0; - self.vy[[i, k]] = 0.5 * (psi.f[[i + 1, k]] + psi.f[[i + 1, k - 1]]); + self.vy[[i, k]] = 0.5 * (psi[[i + 1, k]] + psi[[i + 1, k - 1]]); let i = NX - 1; - self.vy[[i, k]] = -0.5 * (psi.f[[i - 1, k]] + psi.f[[i - 1, k - 1]]); + self.vy[[i, k]] = -0.5 * (psi[[i - 1, k]] + psi[[i - 1, k - 1]]); } } } diff --git a/src/main.rs b/src/main.rs index fff7660..a613a53 100644 --- a/src/main.rs +++ b/src/main.rs @@ -1,5 +1,5 @@ const NY: usize = 31; -const L: f64 = 2.0; +const L: f64 = 2.16; const HEIGHT: f64 = 1.0; const H: f64 = HEIGHT / ((NY - 1) as f64); const NX: usize = ((L / H).round() as usize) + 1; @@ -11,7 +11,7 @@ const PEREODIC: bool = true; use ndarray_stats::QuantileExt; pub mod equations; use csv::Writer; -use equations::{Params, PhiEquation, TemperatureEquation}; +use equations::{ConcentrationEquation, Params, PhiEquation, TemperatureEquation}; use pdifflib::finit_diff::poisson_relax; use serde::Serialize; use std::fs; @@ -27,6 +27,7 @@ struct S { conc: Field2D, temp_eq: TemperatureEquation, phi_eq: PhiEquation, + conc_eq: ConcentrationEquation, params: Params, time: f64, // file:hdf5::File, @@ -35,11 +36,11 @@ impl System for S { fn next_step(&mut self, dt: f64, time: f64) { self.time = time; self.phi_eq - .step(&mut self.phi.f, &self.psi.f, &self.temp.f, dt); + .step(&mut self.phi.f, &self.psi.f, &self.temp.f, &self.conc.f, dt); poisson_relax(&self.phi.f, &mut self.psi.f, H, PEREODIC); self.temp_eq.step(&mut self.temp.f, &self.psi.f, dt); - // self.conc_eq - // .step(&mut self.conc, &self.psi, &self.temp, _dt); + self.conc_eq + .step(&mut self.conc.f, &self.psi.f, &self.temp.f, dt); self.boundary_condition(); } fn fields(&self) -> Vec<&Field2D> { @@ -108,7 +109,7 @@ impl System for S { fn get_max_time(&self) -> f64 { return self.params.time; } - fn get_DT(&self) -> f64 { + fn get_dt(&self) -> f64 { return H * H / 5.0 / 60.; } fn initial_condition(&mut self) { @@ -164,8 +165,9 @@ fn main() { phi: Field2D::new("phi", NX, NY), psi: Field2D::new("psi", NX, NY), conc: Field2D::new("C", NX, NY), - temp_eq: TemperatureEquation::new(NX, NY), - phi_eq: PhiEquation::new(NX, NY, params.rel, params.pr), + temp_eq: TemperatureEquation::new(NX, NY, params.le), + phi_eq: PhiEquation::new(NX, NY, params.rel, params.rel_c, params.sc, params.le), + conc_eq: ConcentrationEquation::new(NX, NY, params.pe), params, time: 0.0, // file:file