mulimoen
/
SummationByParts
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

solve 1D maxwell (periodic)

This commit is contained in:
Magnus Ulimoen 2019-08-09 17:46:59 +02:00
parent 8d51e69607
commit 3ab7e532c6
1 changed files with 121 additions and 73 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

194
src/lib.rs
View File

@ -10,9 +10,9 @@ static ALLOC: wee_alloc::WeeAlloc = wee_alloc::WeeAlloc::INIT;
#[wasm_bindgen] #[wasm_bindgen]
pub struct Universe { pub struct Universe {
e_x: Array2<f32>, ex: Array2<f32>,
// e_y: Array2<f32>, // e_y: Array2<f32>,
// h_z: Array2<f32>, hz: Array2<f32>,
} }
#[wasm_bindgen] #[wasm_bindgen]
@ -36,18 +36,29 @@ fn exponential(x: f32, y: f32, _t: f32) -> f32 {
1.0 / (2.0 * f32::consts::PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp() 1.0 / (2.0 * f32::consts::PI * sigma * sigma) * (-(x * x + y * y) / (2.0 * sigma * sigma)).exp()
} }
fn gaussian(x: f32, y: f32, t: f32) -> f32 {
use std::f32;
let x = x - 0.5;
let y = y - 0.5;
let sigma = 0.05;
1.0 / (2.0 * f32::consts::PI * sigma * sigma)
* (-(x + y - t) * (x + y - t) / (2.0 * sigma * sigma)).exp()
}
#[wasm_bindgen] #[wasm_bindgen]
impl Universe { impl Universe {
pub fn new(width: u32, height: u32) -> Self { pub fn new(width: u32, height: u32) -> Self {
let field = Array2::zeros((height as usize, width as usize)); let field = Array2::zeros((height as usize, width as usize));
let e_x = field.clone(); let ex = field.clone();
// let e_y = field.clone(); // let e_y = field.clone();
// let h_z = field.clone(); let hz = field.clone();
Universe { Universe {
e_x, ex,
// e_y, // e_y,
// h_z, hz,
} }
} }
@ -55,133 +66,162 @@ impl Universe {
let func = match selected_function { let func = match selected_function {
"sin+cos" => sin_plus_cos, "sin+cos" => sin_plus_cos,
"exp" => exponential, "exp" => exponential,
"gauss" => gaussian,
_ => |_x: f32, _y: f32, _t: f32| 0.0, _ => |_x: f32, _y: f32, _t: f32| 0.0,
}; };
let nx = self.e_x.shape()[1]; let nx = self.ex.shape()[1];
let ny = self.e_x.shape()[0]; let ny = self.ex.shape()[0];
for j in 0..ny { for j in 0..ny {
for i in 0..nx { for i in 0..nx {
// Must divice interval on nx/ny instead of nx - 1/ny-1 // Must divice interval on nx/ny instead of nx - 1/ny-1
// due to periodic conditions [0, 1) // due to periodic conditions [0, 1)
let x = i as f32 / nx as f32; let x = i as f32 / nx as f32;
let y = j as f32 / ny as f32; let y = j as f32 / ny as f32;
self.e_x[(j, i)] = func(x, y, t); self.ex[(j, i)] = func(x, y, t);
self.hz[(j, i)] = -func(x, y, t);
// self.e_y[(j as usize, i as usize)] = func(x, y, t); // self.e_y[(j as usize, i as usize)] = func(x, y, t);
// self.h_z[(j as usize, i as usize)] = func(x, y, t);
} }
} }
} }
pub fn advance(&self, fut: &mut Universe, dt: f32, work_buffers: Option<WorkBuffers>) { pub fn advance(&self, fut: &mut Universe, dt: f32, work_buffers: Option<WorkBuffers>) {
assert_eq!(self.e_x.shape(), fut.e_x.shape()); assert_eq!(self.ex.shape(), fut.ex.shape());
let mut y = self.e_x.clone(); let mut y_ex = self.ex.clone();
let mut y_hz = self.hz.clone();
let mut buffers = work_buffers let mut buffers = work_buffers
.unwrap_or_else(|| WorkBuffers::new(self.e_x.shape()[1], self.e_x.shape()[0])); .unwrap_or_else(|| WorkBuffers::new(self.ex.shape()[1], self.ex.shape()[0]));
// y.assign(&self.e_x) let mut k0_ex = &mut buffers.k1_ex;
let mut k0 = &mut buffers.k1; // Only need k0 for first step in RK let mut k0_hz = &mut buffers.k1_hz;
k0.fill(0.0); k0_ex.fill(0.0);
diffx_periodic(y.view(), k0.view_mut()); k0_hz.fill(0.0);
diffy_periodic(y.view(), k0.view_mut()); // u_x
diffx_periodic(y_ex.view(), k0_ex.view_mut());
diffx_periodic(y_hz.view(), k0_hz.view_mut());
// y.assign(&self.e_x); // y = y + 1/2*dt Au_x
y.scaled_add(-1.0 / 2.0 * dt, &k0); y_ex.scaled_add(1.0 / 2.0 * dt, &k0_hz);
let mut k1 = buffers.k1; y_hz.scaled_add(1.0 / 2.0 * dt, &k0_ex);
k1.fill(0.0); let mut k1_ex = buffers.k1_ex;
diffx_periodic(y.view(), k1.view_mut()); let mut k1_hz = buffers.k1_hz;
diffy_periodic(y.view(), k1.view_mut()); k1_ex.fill(0.0);
k1_hz.fill(0.0);
diffx_periodic(y_ex.view(), k1_ex.view_mut());
diffx_periodic(y_hz.view(), k1_hz.view_mut());
y.assign(&self.e_x); y_ex.assign(&self.ex);
y.scaled_add(-1.0 / 2.0 * dt, &k1); y_hz.assign(&self.hz);
let mut k2 = buffers.k2; y_ex.scaled_add(1.0 / 2.0 * dt, &k1_hz);
k2.fill(0.0); y_hz.scaled_add(1.0 / 2.0 * dt, &k1_ex);
diffx_periodic(y.view(), k2.view_mut()); let mut k2_ex = buffers.k2_ex;
diffy_periodic(y.view(), k2.view_mut()); let mut k2_hz = buffers.k2_hz;
k2_ex.fill(0.0);
k2_hz.fill(0.0);
diffx_periodic(y_ex.view(), k2_ex.view_mut());
diffx_periodic(y_hz.view(), k2_hz.view_mut());
y.assign(&self.e_x); // y = y0
y.scaled_add(-1.0 / 2.0 * dt, &k2); y_ex.assign(&self.ex);
let mut k3 = buffers.k3; y_hz.assign(&self.hz);
k3.fill(0.0); y_ex.scaled_add(1.0 / 2.0 * dt, &k2_hz);
diffx_periodic(y.view(), k3.view_mut()); y_hz.scaled_add(1.0 / 2.0 * dt, &k2_ex);
diffy_periodic(y.view(), k3.view_mut()); let mut k3_ex = buffers.k3_ex;
let mut k3_hz = buffers.k3_hz;
k3_ex.fill(0.0);
k3_hz.fill(0.0);
diffx_periodic(y_ex.view(), k3_ex.view_mut());
diffx_periodic(y_hz.view(), k3_hz.view_mut());
y.assign(&self.e_x); y_ex.assign(&self.ex);
y.scaled_add(-dt, &k2); y_hz.assign(&self.hz);
let mut k4 = buffers.k4; y_ex.scaled_add(dt, &k2_hz);
k4.fill(0.0); y_hz.scaled_add(dt, &k2_ex);
diffx_periodic(y.view(), k4.view_mut()); let mut k4_ex = buffers.k4_ex;
diffy_periodic(y.view(), k4.view_mut()); let mut k4_hz = buffers.k4_hz;
k4_ex.fill(0.0);
k4_hz.fill(0.0);
diffx_periodic(y_ex.view(), k4_ex.view_mut());
diffx_periodic(y_hz.view(), k4_hz.view_mut());
Zip::from(&mut fut.e_x) Zip::from(&mut fut.ex)
.and(&self.e_x) .and(&self.ex)
.and(&k1) .and(&k1_hz)
.and(&k2) .and(&k2_hz)
.and(&k3) .and(&k3_hz)
.and(&k4) .and(&k4_hz)
.apply(|y1, &y0, &k1, &k2, &k3, &k4| { .apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + -dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4) *y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
});
Zip::from(&mut fut.hz)
.and(&self.hz)
.and(&k1_ex)
.and(&k2_ex)
.and(&k3_ex)
.and(&k4_ex)
.apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
}); });
} }
pub fn advance_upwind(&self, fut: &mut Universe, dt: f32) { pub fn advance_upwind(&self, fut: &mut Universe, dt: f32) {
assert_eq!(self.e_x.dim(), fut.e_x.dim()); assert_eq!(self.ex.dim(), fut.ex.dim());
let mut ad = Array2::zeros(fut.e_x.dim()); let mut ad = Array2::zeros(fut.ex.dim());
let mut y = self.e_x.clone(); let mut y = self.ex.clone();
// y.assign(&self.e_x) // y.assign(&self.ex)
// ad.fill(0.0); // ad.fill(0.0);
let mut k0 = Array2::zeros(self.e_x.dim()); let mut k0 = Array2::zeros(self.ex.dim());
diffx(y.view(), k0.view_mut()); diffx(y.view(), k0.view_mut());
diffy(y.view(), k0.view_mut()); diffy(y.view(), k0.view_mut());
dissx(y.view(), ad.view_mut()); dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut()); dissy(y.view(), ad.view_mut());
k0.scaled_add(-0.5, &ad); k0.scaled_add(-0.5, &ad);
// y.assign(&self.e_x); // y.assign(&self.ex);
ad.fill(0.0); ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k0); y.scaled_add(-1.0 / 2.0 * dt, &k0);
let mut k1 = Array2::zeros(self.e_x.dim()); let mut k1 = Array2::zeros(self.ex.dim());
diffx(y.view(), k1.view_mut()); diffx(y.view(), k1.view_mut());
diffy(y.view(), k1.view_mut()); diffy(y.view(), k1.view_mut());
dissx(y.view(), ad.view_mut()); dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut()); dissy(y.view(), ad.view_mut());
k1.scaled_add(-0.5, &ad); k1.scaled_add(-0.5, &ad);
y.assign(&self.e_x); y.assign(&self.ex);
ad.fill(0.0); ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k1); y.scaled_add(-1.0 / 2.0 * dt, &k1);
let mut k2 = Array2::zeros(self.e_x.dim()); let mut k2 = Array2::zeros(self.ex.dim());
diffx(y.view(), k2.view_mut()); diffx(y.view(), k2.view_mut());
diffy(y.view(), k2.view_mut()); diffy(y.view(), k2.view_mut());
dissx(y.view(), ad.view_mut()); dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut()); dissy(y.view(), ad.view_mut());
k2.scaled_add(-0.5, &ad); k2.scaled_add(-0.5, &ad);
y.assign(&self.e_x); y.assign(&self.ex);
ad.fill(0.0); ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k2); y.scaled_add(-1.0 / 2.0 * dt, &k2);
let mut k3 = Array2::zeros(self.e_x.dim()); let mut k3 = Array2::zeros(self.ex.dim());
diffx(y.view(), k3.view_mut()); diffx(y.view(), k3.view_mut());
diffy(y.view(), k3.view_mut()); diffy(y.view(), k3.view_mut());
dissx(y.view(), ad.view_mut()); dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut()); dissy(y.view(), ad.view_mut());
k3.scaled_add(-0.5, &ad); k3.scaled_add(-0.5, &ad);
y.assign(&self.e_x); y.assign(&self.ex);
ad.fill(0.0); ad.fill(0.0);
y.scaled_add(-dt, &k2); y.scaled_add(-dt, &k2);
let mut k4 = Array2::zeros(self.e_x.dim()); let mut k4 = Array2::zeros(self.ex.dim());
diffx(y.view(), k4.view_mut()); diffx(y.view(), k4.view_mut());
diffy(y.view(), k4.view_mut()); diffy(y.view(), k4.view_mut());
dissx(y.view(), ad.view_mut()); dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut()); dissy(y.view(), ad.view_mut());
k4.scaled_add(-0.5, &ad); k4.scaled_add(-0.5, &ad);
Zip::from(&mut fut.e_x) Zip::from(&mut fut.ex)
.and(&self.e_x) .and(&self.ex)
.and(&k1) .and(&k1)
.and(&k2) .and(&k2)
.and(&k3) .and(&k3)
@ -192,26 +232,34 @@ impl Universe {
} }
pub fn get_ptr(&mut self) -> *mut u8 { pub fn get_ptr(&mut self) -> *mut u8 {
self.e_x.as_mut_ptr() as *mut u8 self.ex.as_mut_ptr() as *mut u8
} }
} }
#[wasm_bindgen] #[wasm_bindgen]
pub struct WorkBuffers { pub struct WorkBuffers {
k1: Array2<f32>, k1_ex: Array2<f32>,
k2: Array2<f32>, k1_hz: Array2<f32>,
k3: Array2<f32>, k2_ex: Array2<f32>,
k4: Array2<f32>, k2_hz: Array2<f32>,
k3_ex: Array2<f32>,
k3_hz: Array2<f32>,
k4_ex: Array2<f32>,
k4_hz: Array2<f32>,
} }
#[wasm_bindgen] #[wasm_bindgen]
impl WorkBuffers { impl WorkBuffers {
pub fn new(nx: usize, ny: usize) -> Self { pub fn new(nx: usize, ny: usize) -> Self {
Self { Self {
k1: Array2::zeros((ny, nx)), k1_ex: Array2::zeros((ny, nx)),
k2: Array2::zeros((ny, nx)), k1_hz: Array2::zeros((ny, nx)),
k3: Array2::zeros((ny, nx)), k2_ex: Array2::zeros((ny, nx)),
k4: Array2::zeros((ny, nx)), k2_hz: Array2::zeros((ny, nx)),
k3_ex: Array2::zeros((ny, nx)),
k3_hz: Array2::zeros((ny, nx)),
k4_ex: Array2::zeros((ny, nx)),
k4_hz: Array2::zeros((ny, nx)),
} }
} }
} }