firts commit

2026-07-12 15:57:34 +05:00
commit 6e7e8c7815
7 changed files with 1158 additions and 0 deletions
@@ -0,0 +1,184 @@
+use crate::{H, NX, NY, PEREODIC};
+use ndarray::Array2;
+use pdifflib::field::{Field2D, 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};
+
+#[derive(Serialize, Deserialize, Copy, Clone, Debug)]
+pub struct Params {
+    pub rel: f64,
+    pub rel_c: f64,
+    pub le: f64,
+    pub pr: f64,
+    pub pe: f64,
+    pub ma: f64,
+    pub time: f64,
+}
+
+pub struct TemperatureEquation {
+    grid: Grid,
+    delta: Array2<f64>,
+}
+
+impl TemperatureEquation {
+    pub fn new(nx: usize, ny: usize) -> Self {
+        Self {
+            grid: Grid::new(nx, ny),
+            delta: Array2::<f64>::zeros((nx, ny)),
+        }
+    }
+
+    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));
+        }
+        *temp += &self.delta;
+    }
+}
+
+pub struct PhiEquation {
+    grid: Grid,
+    delta: Array2<f64>,
+    rel: f64,
+    pr: f64,
+}
+
+impl PhiEquation {
+    pub fn new(nx: usize, ny: usize, rel: f64, pr: f64) -> Self {
+        Self {
+            grid: Grid::new(nx, ny),
+            delta: Array2::<f64>::zeros((nx, ny)),
+            rel,
+            pr,
+        }
+    }
+
+    pub fn step(
+        &mut self,
+        phi: &mut Array2<f64>,
+        psi: &Array2<f64>,
+        temp: &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 tmp = pr * (lap!(phi, i, j) * inv_h2 + archim * inv_h)
+                + jacobian!(psi, phi, i, j) * inv_h2;
+            self.delta[[i, j]] = tmp * dt;
+        }
+        *phi += &self.delta;
+    }
+}
+
+pub struct PsiEquation;
+
+impl PsiEquation {
+    pub fn new() -> Self {
+        Self
+    }
+
+    pub fn step(&mut self, psi: &mut Array2<f64>, phi: &Array2<f64>) {
+        poisson_relax(&phi, psi, H, PEREODIC);
+    }
+}
+
+pub struct ConcentrationEquation {
+    grid: Grid,
+    qew: Array2<f64>,
+    qsn: Array2<f64>,
+    vx: Array2<f64>,
+    vy: Array2<f64>,
+    pe: f64,
+}
+
+impl ConcentrationEquation {
+    pub fn new(nx: usize, ny: usize, pe: f64) -> Self {
+        Self {
+            grid: Grid::new(nx, ny),
+            qew: Array2::<f64>::zeros((nx, ny)),
+            qsn: Array2::<f64>::zeros((nx, ny)),
+            vx: Array2::<f64>::zeros((nx, ny)),
+            vy: Array2::<f64>::zeros((nx, ny)),
+            pe,
+        }
+    }
+
+    pub fn step(&mut self, conc: &mut Field2D, psi: &Field2D, _temp: &Field2D, 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);
+            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);
+            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;
+        }
+        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]]);
+
+            let k = NY - 1;
+            conc.f[[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]];
+            }
+        } 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]]);
+
+                let i = NX - 1;
+                conc.f[[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]] +=
+                dt / H * 2.0 * (self.qew[[NX - 1, NY - 1]] + self.qsn[[NX - 1, NY - 1]]);
+        }
+    }
+
+    fn calc_v(&mut self, psi: &Field2D) {
+        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]]);
+
+            for k in 1..NY - 1 {
+                self.vx[[i, k]] = 0.5 * (dy!(psi.f, i, k) + dy!(psi.f, i - 1, k));
+            }
+            k = NY - 1;
+            self.vx[[i, k]] = -0.5 * (psi.f[[i, k - 1]] + psi.f[[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));
+            }
+        }
+
+        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]]);
+
+            let i = NX - 1;
+            self.vy[[i, k]] = -0.5 * (psi.f[[i - 1, k]] + psi.f[[i - 1, k - 1]]);
+        }
+    }
+}