From 7a2987126628d74f1ba715f0e288392b57e3b739 Mon Sep 17 00:00:00 2001 From: Magnus Ulimoen Date: Sun, 21 Jul 2019 21:13:54 +0200 Subject: [PATCH] add upwind4 --- src/lib.rs | 169 ++++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 167 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index f67cd91..7de3fab 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -101,13 +101,24 @@ impl Universe { fn diffx(prev: ArrayView2, mut fut: ArrayViewMut2) { for j in 0..prev.shape()[0] { - trad4(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..))); + upwind4(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..))); } } fn diffy(prev: ArrayView2, mut fut: ArrayViewMut2) { for i in 0..prev.shape()[1] { - trad4(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i))); + upwind4(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i))); + } +} + +fn dissx(prev: ArrayView2, mut fut: ArrayViewMut2) { + for j in 0..prev.shape()[0] { + upwind4_diss(prev.slice(s!(j, ..)), fut.slice_mut(s!(j, ..))); + } +} +fn dissy(prev: ArrayView2, mut fut: ArrayViewMut2) { + for i in 0..prev.shape()[1] { + upwind4_diss(prev.slice(s!(.., i)), fut.slice_mut(s!(.., i))); } } @@ -152,3 +163,157 @@ fn trad4(prev: ArrayView1, mut fut: ArrayViewMut1) { + diag[4] * prev[(1)]; fut[(nx - 1)] += diff / dx; } + +fn upwind4(prev: ArrayView1, mut fut: ArrayViewMut1) { + assert_eq!(prev.shape(), fut.shape()); + let nx = prev.shape()[0]; + + let dx = 1.0 / (nx - 1) as f32; + + let diag = [ + -1.0 / 24.0, + 1.0 / 4.0, + -7.0 / 8.0, + 0.0, + 7.0 / 8.0, + -1.0 / 4.0, + 1.0 / 24.0, + ]; + + let diff = diag[0] * prev[(nx - 3)] + + diag[1] * prev[(nx - 2)] + + diag[2] * prev[(nx - 1)] + + diag[3] * prev[(0)] + + diag[4] * prev[(1)] + + diag[5] * prev[(2)] + + diag[6] * prev[(3)]; + fut[0] += diff / dx; + let diff = diag[0] * prev[(nx - 2)] + + diag[1] * prev[(nx - 1)] + + diag[2] * prev[(0)] + + diag[3] * prev[(1)] + + diag[4] * prev[(2)] + + diag[5] * prev[(3)] + + diag[6] * prev[(4)]; + fut[1] += diff / dx; + let diff = diag[0] * prev[(nx - 1)] + + diag[1] * prev[(0)] + + diag[2] * prev[(1)] + + diag[3] * prev[(2)] + + diag[4] * prev[(3)] + + diag[5] * prev[(4)] + + diag[6] * prev[(5)]; + fut[2] += diff / dx; + + for i in 3..nx - 3 { + let diff = diag[0] * prev[(i - 3)] + + diag[1] * prev[(i - 2)] + + diag[2] * prev[(i - 1)] + + diag[3] * prev[(i)] + + diag[4] * prev[(i + 1)] + + diag[5] * prev[(i + 2)] + + diag[6] * prev[(i + 3)]; + fut[(i)] += diff / dx; + } + let diff = diag[0] * prev[(nx - 6)] + + diag[1] * prev[(nx - 5)] + + diag[2] * prev[(nx - 4)] + + diag[3] * prev[(nx - 3)] + + diag[4] * prev[(nx - 2)] + + diag[5] * prev[(nx - 1)] + + diag[6] * prev[(0)]; + fut[(nx - 3)] += diff / dx; + let diff = diag[0] * prev[(nx - 5)] + + diag[1] * prev[(nx - 4)] + + diag[2] * prev[(nx - 3)] + + diag[3] * prev[(nx - 2)] + + diag[4] * prev[(nx - 1)] + + diag[5] * prev[(0)] + + diag[6] * prev[(1)]; + fut[(nx - 2)] += diff / dx; + let diff = diag[0] * prev[(nx - 4)] + + diag[1] * prev[(nx - 3)] + + diag[2] * prev[(nx - 2)] + + diag[3] * prev[(nx - 1)] + + diag[4] * prev[(0)] + + diag[5] * prev[(1)] + + diag[6] * prev[(2)]; + fut[(nx - 1)] += diff / dx; +} + +fn upwind4_diss(prev: ArrayView1, mut fut: ArrayViewMut1) { + assert_eq!(prev.shape(), fut.shape()); + let nx = prev.shape()[0]; + + let dx = 1.0 / (nx - 1) as f32; + + let diag = [ + 1.0 / 24.0, + -1.0 / 4.0, + 5.0 / 8.0, + -5.0 / 6.0, + 5.0 / 8.0, + -1.0 / 4.0, + 1.0 / 24.0, + ]; + + let diff = diag[0] * prev[(nx - 3)] + + diag[1] * prev[(nx - 2)] + + diag[2] * prev[(nx - 1)] + + diag[3] * prev[(0)] + + diag[4] * prev[(1)] + + diag[5] * prev[(2)] + + diag[6] * prev[(3)]; + fut[0] += diff / dx; + let diff = diag[0] * prev[(nx - 2)] + + diag[1] * prev[(nx - 1)] + + diag[2] * prev[(0)] + + diag[3] * prev[(1)] + + diag[4] * prev[(2)] + + diag[5] * prev[(3)] + + diag[6] * prev[(4)]; + fut[1] += diff / dx; + let diff = diag[0] * prev[(nx - 1)] + + diag[1] * prev[(0)] + + diag[2] * prev[(1)] + + diag[3] * prev[(2)] + + diag[4] * prev[(3)] + + diag[5] * prev[(4)] + + diag[6] * prev[(5)]; + fut[2] += diff / dx; + + for i in 3..nx - 3 { + let diff = diag[0] * prev[(i - 3)] + + diag[1] * prev[(i - 2)] + + diag[2] * prev[(i - 1)] + + diag[3] * prev[(i)] + + diag[4] * prev[(i + 1)] + + diag[5] * prev[(i + 2)] + + diag[6] * prev[(i + 3)]; + fut[(i)] += diff / dx; + } + let diff = diag[0] * prev[(nx - 6)] + + diag[1] * prev[(nx - 5)] + + diag[2] * prev[(nx - 4)] + + diag[3] * prev[(nx - 3)] + + diag[4] * prev[(nx - 2)] + + diag[5] * prev[(nx - 1)] + + diag[6] * prev[(0)]; + fut[(nx - 3)] += diff / dx; + let diff = diag[0] * prev[(nx - 5)] + + diag[1] * prev[(nx - 4)] + + diag[2] * prev[(nx - 3)] + + diag[3] * prev[(nx - 2)] + + diag[4] * prev[(nx - 1)] + + diag[5] * prev[(0)] + + diag[6] * prev[(1)]; + fut[(nx - 2)] += diff / dx; + let diff = diag[0] * prev[(nx - 4)] + + diag[1] * prev[(nx - 3)] + + diag[2] * prev[(nx - 2)] + + diag[3] * prev[(nx - 1)] + + diag[4] * prev[(0)] + + diag[5] * prev[(1)] + + diag[6] * prev[(2)]; + fut[(nx - 1)] += diff / dx; +}