fix

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"

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"

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<f64>,
+    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::<f64>::zeros((nx, ny)),
+            le,
         }
     }
 
     pub fn step(&mut self, temp: &mut Array2<f64>, psi: &Array2<f64>, 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<f64>,
     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::<f64>::zeros((nx, ny)),
             rel,
-            pr,
+            rel_c,
+            sc,
+            le,
         }
     }
 
@@ -60,14 +67,14 @@ impl PhiEquation {
         phi: &mut Array2<f64>,
         psi: &Array2<f64>,
         temp: &Array2<f64>,
+        conc: &Array2<f64>,
         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 archim = self.le * self.rel * dx!(temp, i, j) + self.rel_c * dx!(conc, i, j);
-            let tmp = pr * (lap!(phi, i, j) * inv_h2 + archim * inv_h)
+            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<f64>,
+        psi: &Array2<f64>,
+        _temp: &Array2<f64>,
+        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);
+            let mut tmp = self.vx[[i, k]] * mx_b!(conc, i, k);
-            tmp += -dx_b!(conc.f, 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);
+            let mut tmp = ((-self.vy[[i, k]]) + self.pe * H) * my_b!(conc, i, k);
-            tmp += -dy_b!(conc.f, 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[[0, k]] = conc[[NX - 2, k]];
-                conc.f[[NX - 1, k]] = conc.f[[1, 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[[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[[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[[NX - 1, 0]] += dt / H * 2.0 * (self.qew[[NX - 1, 0]] - self.qsn[[NX - 1, 1]]);
-            conc.f[[NX - 1, NY - 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<f64>) {
         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]]);
         }
     }
 }

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
+        self.conc_eq
-        //      .step(&mut self.conc, &self.psi, &self.temp, _dt);
+            .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),
+        temp_eq: TemperatureEquation::new(NX, NY, params.le),
-        phi_eq: PhiEquation::new(NX, NY, params.rel, params.pr),
+        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