From b80ca3f5e8de92b75d5f791a17fab98e6ec22032 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Fri, 18 Sep 2020 17:57:57 +0200 Subject: [PATCH] Add D2 operator and heat equation --- Cargo.toml | 1 + heat-equation/Cargo.toml | 11 +++++ heat-equation/src/main.rs | 68 +++++++++++++++++++++++++++++++ sbp/src/operators.rs | 6 +++ sbp/src/operators/traditional4.rs | 38 +++++++++++++++-- 5 files changed, 120 insertions(+), 4 deletions(-) create mode 100644 heat-equation/Cargo.toml create mode 100644 heat-equation/src/main.rs diff --git a/Cargo.toml b/Cargo.toml index dd18645..12431b2 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,6 +7,7 @@ members = [ "maxwell", "shallow_water", "gridgeneration", + "heat-equation", ] default-members = ["sbp", "euler", "maxwell", "shallow_water"] diff --git a/heat-equation/Cargo.toml b/heat-equation/Cargo.toml new file mode 100644 index 0000000..bc8ecb8 --- /dev/null +++ b/heat-equation/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "heat-equation" +version = "0.1.0" +authors = ["Magnus Ulimoen "] +edition = "2018" + + +[dependencies] +sbp = { path = "../sbp" } +ndarray = "0.13.1" +plotters = { version = "0.3.0", default-features = false, features = ["bitmap_gif", "bitmap_backend", "line_series"] } diff --git a/heat-equation/src/main.rs b/heat-equation/src/main.rs new file mode 100644 index 0000000..cd8d22d --- /dev/null +++ b/heat-equation/src/main.rs @@ -0,0 +1,68 @@ +use ndarray::Array1; +use plotters::prelude::*; +use sbp::{ + integrate::{integrate, Rk4}, + operators::{SbpOperator1d2, SBP4}, + Float, +}; + +fn main() { + let drawing_area = BitMapBackend::gif("result.gif", (300, 300), 100) + .unwrap() + .into_drawing_area(); + let mut chart = ChartBuilder::on(&drawing_area) + .x_label_area_size(40) + .y_label_area_size(40) + .build_cartesian_2d(0.0..1.0, -1.0..2.0) + .unwrap(); + + let nx: usize = 101; + let dt = 0.2 / nx.pow(2) as Float / 3.0; + let x = Array1::from_shape_fn((nx,), |i| i as Float / (nx - 1) as Float); + let v0 = x.map(|&x| 0.0 * x * x * x * x); + + let op = SBP4; + + let mut k = [v0.clone(), v0.clone(), v0.clone(), v0.clone()]; + let rhs = move |fut: &mut Array1, prev: &Array1, _t: Float| { + fut.fill(0.0); + op.diff2(prev.view(), fut.view_mut()); + + let h = 1.0 / (nx - 1) as Float; + let tau = (1.0, -1.0); + let v0 = 1.0; + let vn = -1.0; + + let d1 = op.d1(); + + for (d, fut) in d1.iter().zip(fut.iter_mut()) { + *fut += tau.0 / (h * h) * d * (prev[0] - v0); + } + for (d, fut) in d1.iter().zip(fut.iter_mut().rev().take(d1.len()).rev()) { + *fut += tau.1 / (h * h) * d * (prev[nx - 1] - vn); + } + }; + + let mut v1 = v0.clone(); + let mut v2 = v0; + for i in 0..90 { + if i % 3 == 0 { + drawing_area.fill(&WHITE).unwrap(); + chart + .configure_mesh() + .x_desc("x") + .y_desc("y") + .draw() + .unwrap(); + chart + .draw_series(LineSeries::new( + x.iter().zip(v1.iter()).map(|(&x, &y)| (x, y)), + &BLACK, + )) + .unwrap(); + drawing_area.present().unwrap(); + } + integrate::(rhs, &v1, &mut v2, &mut 0.0, dt, &mut k); + std::mem::swap(&mut v1, &mut v2); + } +} diff --git a/sbp/src/operators.rs b/sbp/src/operators.rs index 180fcf3..67e20a8 100644 --- a/sbp/src/operators.rs +++ b/sbp/src/operators.rs @@ -17,6 +17,12 @@ pub trait SbpOperator1d: Send + Sync { fn h_matrix(&self, n: usize) -> sprs::CsMat; } +pub trait SbpOperator1d2: SbpOperator1d { + fn diff2(&self, prev: ArrayView1, fut: ArrayViewMut1); + /// Lacks a scaling of 1/h^2 + fn d1(&self) -> &[Float]; +} + pub trait SbpOperator2d: Send + Sync { fn diffxi(&self, prev: ArrayView2, fut: ArrayViewMut2); fn diffeta(&self, prev: ArrayView2, fut: ArrayViewMut2); diff --git a/sbp/src/operators/traditional4.rs b/sbp/src/operators/traditional4.rs index f17d5cd..f7495cd 100644 --- a/sbp/src/operators/traditional4.rs +++ b/sbp/src/operators/traditional4.rs @@ -16,10 +16,22 @@ impl SBP4 { ]; #[rustfmt::skip] const BLOCK: &'static [&'static [Float]] = &[ - &[-1.41176470588235e+00, 1.73529411764706e+00, -2.35294117647059e-01, -8.82352941176471e-02], - &[-5.00000000000000e-01, 0.00000000000000e+00, 5.00000000000000e-01], - &[9.30232558139535e-02, -6.86046511627907e-01, 0.00000000000000e+00, 6.86046511627907e-01, -9.30232558139535e-02], - &[3.06122448979592e-02, 0.00000000000000e+00, -6.02040816326531e-01, 0.00000000000000e+00, 6.53061224489796e-01, -8.16326530612245e-02], + &[-24.0/17.0, 59.0/34.0, -4.0/17.0, -3.0/34.0], + &[-1.0/2.0, 0.0, 1.0/2.0], + &[4.0/43.0, -59.0/86.0, 0.0, 59.0/86.0, -4.0/43.0], + &[3.0/98.0, 0.0, -59.0/98.0, 0.0, 32.0/49.0, -4.0/49.0] + ]; + + #[rustfmt::skip] + const D2DIAG: &'static [Float] = &[ + -1.0 / 12.0, 4.0 / 3.0, -5.0 / 2.0, 4.0 / 3.0, -1.0 / 12.0 + ]; + #[rustfmt::skip] + const D2BLOCK: &'static [&'static [Float]] = &[ + &[2.0, -5.0, 4.0, -1.0], + &[1.0, -2.0, 1.0], + &[-4.0/43.0, 59.0/43.0, -110.0/43.0, 59.0/43.0, -4.0/43.0], + &[-1.0/49.0, 0.0, 59.0/49.0, -118.0/49.0, 64.0/49.0, -4.0/49.0] ]; } @@ -81,6 +93,24 @@ impl SbpOperator2d for (&SBP, &SBP4) { } } +impl super::SbpOperator1d2 for SBP4 { + fn diff2(&self, prev: ArrayView1, mut fut: ArrayViewMut1) { + super::diff_op_1d( + Self::D2BLOCK, + Self::D2DIAG, + super::Symmetry::Symmetric, + super::OperatorType::Normal, + prev, + fut.view_mut(), + ); + let hi = (prev.len() - 1) as Float; + fut.map_inplace(|x| *x *= hi) + } + fn d1(&self) -> &[Float] { + &[-88.0 / 17.0, 144.0 / 17.0, -72.0 / 17.0, 16.0 / 17.0] + } +} + #[test] fn test_trad4() { use super::testing::*;