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

add upwind dissipation advance()

This commit is contained in:
Magnus Ulimoen 2019-07-22 17:53:08 +02:00
parent 7a29871266
commit d52c461d6d
2 changed files with 97 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
main.js
View File

@ -90,7 +90,7 @@ async function run() {
let universes = [Universe.new(width, height), Universe.new(width, height)]; let universes = [Universe.new(width, height), Universe.new(width, height)];
let t = performance.now()/1000.0; let t = performance.now()/1000.0;
universes[0].set_initial(t); universes[0].set_initial(t, "sin+cos");
const field = new Float32Array(wasm.memory.buffer, const field = new Float32Array(wasm.memory.buffer,

105
src/lib.rs
View File

@ -1,4 +1,4 @@
use ndarray::{s, Array2, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2}; use ndarray::{s, Array2, ArrayView1, ArrayView2, ArrayViewMut1, ArrayViewMut2, Zip};
use wasm_bindgen::prelude::*; use wasm_bindgen::prelude::*;
#[cfg(feature = "wee_alloc")] #[cfg(feature = "wee_alloc")]
@ -20,11 +20,21 @@ pub fn set_panic_hook() {
console_error_panic_hook::set_once(); console_error_panic_hook::set_once();
} }
fn func(x: f32, y: f32, t: f32) -> f32 { fn sin_plus_cos(x: f32, y: f32, t: f32) -> f32 {
use std::f32; use std::f32;
(2.0 * f32::consts::PI * (x + y) - t).sin() (2.0 * f32::consts::PI * (x + y) - t).sin()
} }
fn exponential(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 * x + y * y) / (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 {
@ -42,7 +52,12 @@ impl Universe {
} }
} }
pub fn set_initial(&mut self, t: f32) { pub fn set_initial(&mut self, t: f32, selected_function: &str) {
let func = match selected_function {
"sin+cos" => sin_plus_cos,
"exp" => exponential,
_ => |_x: f32, _y: f32, _t: f32| 0.0,
};
for j in 0..self.height { for j in 0..self.height {
for i in 0..self.width { for i in 0..self.width {
let x = i as f32 / self.width as f32; let x = i as f32 / self.width as f32;
@ -53,12 +68,13 @@ impl Universe {
} }
} }
} }
pub fn advance(&self, fut: &mut Universe, dt: f32) { pub fn advance(&self, fut: &mut Universe, dt: f32) {
assert_eq!(self.width, fut.width); assert_eq!(self.width, fut.width);
assert_eq!(self.height, fut.height); assert_eq!(self.height, fut.height);
let mut y = self.e_x.clone(); let mut y = self.e_x.clone();
// y.assign(&self.e_x)
let mut k0 = Array2::zeros((self.height as usize, self.width as usize)); let mut k0 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k0.view_mut()); diffx(y.view(), k0.view_mut());
diffy(y.view(), k0.view_mut()); diffy(y.view(), k0.view_mut());
@ -87,11 +103,82 @@ impl Universe {
diffx(y.view(), k4.view_mut()); diffx(y.view(), k4.view_mut());
diffy(y.view(), k4.view_mut()); diffy(y.view(), k4.view_mut());
fut.e_x.assign(&self.e_x); Zip::from(&mut fut.e_x)
fut.e_x.scaled_add(-dt / 6.0, &k1); .and(&self.e_x)
fut.e_x.scaled_add(-2.0 * dt / 6.0, &k2); .and(&k1)
fut.e_x.scaled_add(-2.0 * dt / 6.0, &k3); .and(&k2)
fut.e_x.scaled_add(-dt / 6.0, &k4); .and(&k3)
.and(&k4)
.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) {
assert_eq!(self.width, fut.width);
assert_eq!(self.height, fut.height);
let mut ad = Array2::zeros(fut.e_x.dim());
let mut y = self.e_x.clone();
// y.assign(&self.e_x)
// ad.fill(0.0);
let mut k0 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k0.view_mut());
diffy(y.view(), k0.view_mut());
dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut());
k0.scaled_add(-0.5, &ad);
// y.assign(&self.e_x);
ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k0);
let mut k1 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k1.view_mut());
diffy(y.view(), k1.view_mut());
dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut());
k1.scaled_add(-0.5, &ad);
y.assign(&self.e_x);
ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k1);
let mut k2 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k2.view_mut());
diffy(y.view(), k2.view_mut());
dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut());
k2.scaled_add(-0.5, &ad);
y.assign(&self.e_x);
ad.fill(0.0);
y.scaled_add(-1.0 / 2.0 * dt, &k2);
let mut k3 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k3.view_mut());
diffy(y.view(), k3.view_mut());
dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut());
k3.scaled_add(-0.5, &ad);
y.assign(&self.e_x);
ad.fill(0.0);
y.scaled_add(-dt, &k2);
let mut k4 = Array2::zeros((self.height as usize, self.width as usize));
diffx(y.view(), k4.view_mut());
diffy(y.view(), k4.view_mut());
dissx(y.view(), ad.view_mut());
dissy(y.view(), ad.view_mut());
k4.scaled_add(-0.5, &ad);
Zip::from(&mut fut.e_x)
.and(&self.e_x)
.and(&k1)
.and(&k2)
.and(&k3)
.and(&k4)
.apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + -dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
});
} }
pub fn get_ptr(&mut self) -> *mut u8 { pub fn get_ptr(&mut self) -> *mut u8 {