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

solve 2D maxwell

This commit is contained in:
Magnus Ulimoen 2019-08-09 21:00:15 +02:00
parent 6816d6bd43
commit 43f1088d66
2 changed files with 73 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

50
main.js
View File

@ -33,6 +33,7 @@ async function run() {
varying lowp vec2 vVertexPosition; varying lowp vec2 vVertexPosition;
uniform sampler2D uSamplerEX; uniform sampler2D uSamplerEX;
uniform sampler2D uSamplerEY;
uniform sampler2D uSamplerHZ; uniform sampler2D uSamplerHZ;
uniform int uChosenField; uniform int uChosenField;
@ -43,6 +44,9 @@ async function run() {
if (uChosenField == 1) { if (uChosenField == 1) {
r = texture2D(uSamplerEX, vVertexPosition).a; r = texture2D(uSamplerEX, vVertexPosition).a;
r = (r + 1.0)/2.0; r = (r + 1.0)/2.0;
} else if (uChosenField == 3) {
b = texture2D(uSamplerEY, vVertexPosition).a;
b = (b + 1.0)/2.0;
} else if (uChosenField == 2) { } else if (uChosenField == 2) {
g = texture2D(uSamplerHZ, vVertexPosition).a; g = texture2D(uSamplerHZ, vVertexPosition).a;
g = (g + 1.0)/2.0; g = (g + 1.0)/2.0;
@ -99,32 +103,31 @@ async function run() {
gl.enableVertexAttribArray(attrib_vertex_location); gl.enableVertexAttribArray(attrib_vertex_location);
} }
const texture_ex = gl.createTexture(); const create_2D_texture = function(loc, name) {
gl.activeTexture(gl.TEXTURE0); const texture = gl.createTexture();
gl.bindTexture(gl.TEXTURE_2D, texture_ex); gl.activeTexture(gl.TEXTURE0 + loc);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.bindTexture(gl.TEXTURE_2D, texture);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
const uniform_sampler_ex = gl.getUniformLocation(shaderProgram, 'uSamplerEX'); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
gl.uniform1i(uniform_sampler_ex, 0); const uniform_sampler = gl.getUniformLocation(shaderProgram, name);
gl.uniform1i(uniform_sampler, loc);
const texture_hz = gl.createTexture(); return texture;
gl.activeTexture(gl.TEXTURE1); }
gl.bindTexture(gl.TEXTURE_2D, texture_hz); const texture_ex = create_2D_texture(0, 'uSamplerEX');
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); const texture_ey = create_2D_texture(2, 'uSamplerEY');
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); const texture_hz = create_2D_texture(1, 'uSamplerHZ');
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
const uniform_sampler_hz = gl.getUniformLocation(shaderProgram, 'uSamplerHZ');
gl.uniform1i(uniform_sampler_hz, 1);
const chosen_field = { const chosen_field = {
uLocation: gl.getUniformLocation(shaderProgram, 'uChosenField'), uLocation: gl.getUniformLocation(shaderProgram, 'uChosenField'),
value: 2, value: 3,
cycle: function() { cycle: function() {
if (this.value == 1) { if (this.value == 1) {
this.value = 2; this.value = 2;
} else if (this.value == 2) {
this.value = 3;
} else { } else {
this.value = 1; this.value = 1;
} }
@ -134,7 +137,7 @@ async function run() {
chosen_field.cycle(); chosen_field.cycle();
gl.clearColor(1.0, 0.753, 0.796, 1.0); gl.clearColor(1.0, 0.753, 0.796, 1.0); // A nice pink
gl.clearDepth(1.0); gl.clearDepth(1.0);
gl.enable(gl.DEPTH_TEST); gl.enable(gl.DEPTH_TEST);
gl.depthFunc(gl.LEQUAL); gl.depthFunc(gl.LEQUAL);
@ -150,7 +153,7 @@ async function run() {
const TIMEFACTOR = 1.0/7000; const TIMEFACTOR = 1.0/7000;
let t = performance.now()*TIMEFACTOR; let t = performance.now()*TIMEFACTOR;
universes[0].set_initial(t, "sin+cos"); universes[0].set_initial(t, "exp");
function drawMe(t_draw) { function drawMe(t_draw) {
gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT); gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT);
@ -168,6 +171,9 @@ async function run() {
const field_ex = new Float32Array(wasm.memory.buffer, const field_ex = new Float32Array(wasm.memory.buffer,
universes[0].get_ex_ptr(), universes[0].get_ex_ptr(),
width*height); width*height);
const field_ey = new Float32Array(wasm.memory.buffer,
universes[0].get_ey_ptr(),
width*height);
const field_hz = new Float32Array(wasm.memory.buffer, const field_hz = new Float32Array(wasm.memory.buffer,
universes[0].get_hz_ptr(), universes[0].get_hz_ptr(),
width*height); width*height);
@ -179,6 +185,8 @@ async function run() {
const srcType = gl.FLOAT; const srcType = gl.FLOAT;
gl.activeTexture(gl.TEXTURE0); gl.activeTexture(gl.TEXTURE0);
gl.texImage2D(gl.TEXTURE_2D, level, internalFormat, width, height, border, srcFormat, srcType, field_ex); gl.texImage2D(gl.TEXTURE_2D, level, internalFormat, width, height, border, srcFormat, srcType, field_ex);
gl.activeTexture(gl.TEXTURE2);
gl.texImage2D(gl.TEXTURE_2D, level, internalFormat, width, height, border, srcFormat, srcType, field_ey);
gl.activeTexture(gl.TEXTURE1); gl.activeTexture(gl.TEXTURE1);
gl.texImage2D(gl.TEXTURE_2D, level, internalFormat, width, height, border, srcFormat, srcType, field_hz); gl.texImage2D(gl.TEXTURE_2D, level, internalFormat, width, height, border, srcFormat, srcType, field_hz);
} }

78
src/lib.rs
View File

@ -2,7 +2,7 @@ use ndarray::{Array2, Zip};
use wasm_bindgen::prelude::*; use wasm_bindgen::prelude::*;
mod operators; mod operators;
use operators::diffx_periodic; use operators::{diffx_periodic, diffy_periodic};
#[cfg(feature = "wee_alloc")] #[cfg(feature = "wee_alloc")]
#[global_allocator] #[global_allocator]
@ -11,7 +11,7 @@ static ALLOC: wee_alloc::WeeAlloc = wee_alloc::WeeAlloc::INIT;
#[wasm_bindgen] #[wasm_bindgen]
pub struct Universe { pub struct Universe {
ex: Array2<f32>, ex: Array2<f32>,
// e_y: Array2<f32>, ey: Array2<f32>,
hz: Array2<f32>, hz: Array2<f32>,
} }
@ -36,37 +36,21 @@ 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 ex = field.clone(); let ex = field.clone();
// let e_y = field.clone(); let ey = field.clone();
let hz = field.clone(); let hz = field;
Universe { Universe { ex, ey, hz }
ex,
// e_y,
hz,
}
} }
pub fn set_initial(&mut self, t: f32, selected_function: &str) { pub fn set_initial(&mut self, t: f32, selected_function: &str) {
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.ex.shape()[1]; let nx = self.ex.shape()[1];
@ -77,9 +61,9 @@ impl Universe {
// 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.ex[(j, i)] = func(x, y, t); self.ex[(j, i)] = 0.0;
self.hz[(j, i)] = -func(x, y, t); self.ey[(j, i)] = 0.0;
// self.e_y[(j as usize, i as usize)] = func(x, y, t); self.hz[(j, i)] = func(x, y, t);
} }
} }
} }
@ -97,19 +81,23 @@ impl Universe {
// y = y0 + c*kn // y = y0 + c*kn
y.0.assign(&self.ex); y.0.assign(&self.ex);
y.1.assign(&self.hz); y.1.assign(&self.hz);
y.2.assign(&self.ey);
match i { match i {
0 => {} 0 => {}
1 => { 1 => {
y.0.scaled_add(1.0 / 2.0 * dt, &k[i - 1].0); y.0.scaled_add(1.0 / 2.0 * dt, &k[i - 1].0);
y.1.scaled_add(1.0 / 2.0 * dt, &k[i - 1].1); y.1.scaled_add(1.0 / 2.0 * dt, &k[i - 1].1);
y.2.scaled_add(1.0 / 2.0 * dt, &k[i - 1].2);
} }
2 => { 2 => {
y.0.scaled_add(1.0 / 2.0 * dt, &k[i - 1].0); y.0.scaled_add(1.0 / 2.0 * dt, &k[i - 1].0);
y.1.scaled_add(1.0 / 2.0 * dt, &k[i - 1].1); y.1.scaled_add(1.0 / 2.0 * dt, &k[i - 1].1);
y.2.scaled_add(1.0 / 2.0 * dt, &k[i - 1].2);
} }
3 => { 3 => {
y.0.scaled_add(dt, &k[i - 1].0); y.0.scaled_add(dt, &k[i - 1].0);
y.1.scaled_add(dt, &k[i - 1].1); y.1.scaled_add(dt, &k[i - 1].1);
y.2.scaled_add(dt, &k[i - 1].2);
} }
_ => { _ => {
unreachable!(); unreachable!();
@ -117,10 +105,19 @@ impl Universe {
}; };
k[i].0.fill(0.0); k[i].0.fill(0.0);
k[i].1.fill(0.0); // ex = hz_y
diffy_periodic(y.1.view(), k[i].0.view_mut());
diffx_periodic(y.0.view(), k[i].1.view_mut()); // ey = -hz_x
diffx_periodic(y.1.view(), k[i].0.view_mut()); k[i].2.fill(0.0);
diffx_periodic(y.1.view(), k[i].2.view_mut());
k[i].2.mapv_inplace(|v| -v);
// hz = -ey_x + ex_y
k[i].1.fill(0.0);
diffx_periodic(y.2.view(), k[i].1.view_mut());
k[i].1.mapv_inplace(|v| -v);
diffy_periodic(y.0.view(), k[i].1.view_mut());
} }
Zip::from(&mut fut.ex) Zip::from(&mut fut.ex)
@ -141,12 +138,25 @@ impl Universe {
.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.ey)
.and(&self.ey)
.and(&k[0].2)
.and(&k[1].2)
.and(&k[2].2)
.and(&k[3].2)
.apply(|y1, &y0, &k1, &k2, &k3, &k4| {
*y1 = y0 + dt / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
});
} }
pub fn get_ex_ptr(&mut self) -> *mut u8 { pub fn get_ex_ptr(&mut self) -> *mut u8 {
self.ex.as_mut_ptr() as *mut u8 self.ex.as_mut_ptr() as *mut u8
} }
pub fn get_ey_ptr(&mut self) -> *mut u8 {
self.ey.as_mut_ptr() as *mut u8
}
pub fn get_hz_ptr(&mut self) -> *mut u8 { pub fn get_hz_ptr(&mut self) -> *mut u8 {
self.hz.as_mut_ptr() as *mut u8 self.hz.as_mut_ptr() as *mut u8
} }
@ -154,8 +164,8 @@ impl Universe {
#[wasm_bindgen] #[wasm_bindgen]
pub struct WorkBuffers { pub struct WorkBuffers {
y: (Array2<f32>, Array2<f32>), y: (Array2<f32>, Array2<f32>, Array2<f32>),
buf: [(Array2<f32>, Array2<f32>); 4], buf: [(Array2<f32>, Array2<f32>, Array2<f32>); 4],
} }
#[wasm_bindgen] #[wasm_bindgen]
@ -163,12 +173,12 @@ impl WorkBuffers {
pub fn new(nx: usize, ny: usize) -> Self { pub fn new(nx: usize, ny: usize) -> Self {
let arr = Array2::zeros((ny, nx)); let arr = Array2::zeros((ny, nx));
Self { Self {
y: (arr.clone(), arr.clone()), y: (arr.clone(), arr.clone(), arr.clone()),
buf: [ buf: [
(arr.clone(), arr.clone()), (arr.clone(), arr.clone(), arr.clone()),
(arr.clone(), arr.clone()), (arr.clone(), arr.clone(), arr.clone()),
(arr.clone(), arr.clone()), (arr.clone(), arr.clone(), arr.clone()),
(arr.clone(), arr), (arr.clone(), arr.clone(), arr),
], ],
} }
} }