shallow water equations solver
This commit is contained in:
Magnus Ulimoen
10
shallow_water/Cargo.toml
Normal file
10
shallow_water/Cargo.toml
Normal file
@@ -0,0 +1,10 @@
|
||||
[package]
|
||||
name = "shallow-water"
|
||||
version = "0.1.0"
|
||||
authors = ["Magnus Ulimoen <flymagnus@gmail.com>"]
|
||||
edition = "2018"
|
||||
|
||||
[dependencies]
|
||||
ndarray = "0.13.1"
|
||||
sbp = { path = "../sbp" }
|
||||
log = "0.4.8"
|
||||
6
shallow_water/misc/.gitignore
vendored
Normal file
6
shallow_water/misc/.gitignore
vendored
Normal file
@@ -0,0 +1,6 @@
|
||||
*.aux
|
||||
*.fdb_latexmk
|
||||
*.fls
|
||||
*.log
|
||||
*.pdf
|
||||
*.synctex.gz
|
||||
102
shallow_water/misc/fluxes.py
Executable file
102
shallow_water/misc/fluxes.py
Executable file
@@ -0,0 +1,102 @@
|
||||
#! /usr/bin/env python3
|
||||
import sympy as sp
|
||||
from sympy.utilities.codegen import codegen
|
||||
|
||||
x, y, t = sp.symbols("x,y,t")
|
||||
rho, g = sp.symbols("rho,g", positive=True, real=True)
|
||||
# Ignoring rho, assumed to be constant
|
||||
eta = sp.Symbol("eta", positive=True, real=True)
|
||||
etau, etav = sp.symbols("etau,etav", real=True)
|
||||
|
||||
u = etau / eta
|
||||
v = etav / eta
|
||||
q = sp.Matrix([eta, etau, etav])
|
||||
|
||||
E = sp.Matrix([eta * u, eta * u ** 2 + g * eta ** 2 / 2, eta * u * v])
|
||||
F = sp.Matrix([eta * v, eta * u * v, eta * v ** 2 + g * eta ** 2 / 2])
|
||||
|
||||
q = sp.Matrix([eta, etau, etav])
|
||||
A = E.jacobian(q)
|
||||
sp.pprint(A)
|
||||
B = F.jacobian(q)
|
||||
sp.pprint(B)
|
||||
|
||||
f = sp.symbols("f")
|
||||
coriolis = sp.Matrix([0, -f * v * eta, +f * u * eta])
|
||||
sp.pprint(coriolis)
|
||||
b = sp.symbols("b")
|
||||
frictional = sp.Matrix([0, b * u * eta, b * v * eta])
|
||||
sp.pprint(frictional)
|
||||
|
||||
print("A:")
|
||||
# We can also diagonalise the transpose,
|
||||
# which still gives the same (linearised) system.
|
||||
# This results in a much nicer formulation for A+ and A-
|
||||
A = A.T
|
||||
S, D = A.diagonalize()
|
||||
print("Diagonals:")
|
||||
sp.pprint(D[0, 0])
|
||||
sp.pprint(D[1, 1])
|
||||
sp.pprint(D[2, 2])
|
||||
print("S:")
|
||||
sp.pprint(S)
|
||||
print("SI:")
|
||||
sp.pprint(sp.simplify(S.inv()))
|
||||
|
||||
m = abs(etau / eta) + abs(sp.sqrt(eta ** (3)) * sp.sqrt(g) / eta)
|
||||
|
||||
L = D + sp.Matrix([[m, 0, 0], [0, m, 0], [0, 0, m]])
|
||||
Aplus = sp.simplify(S * L * S.inv()) / 2
|
||||
Aplus = Aplus.T
|
||||
L = D - sp.Matrix([[m, 0, 0], [0, m, 0], [0, 0, m]])
|
||||
Aminus = sp.simplify(S * L * S.inv()) / 2
|
||||
Aminus = Aminus.T
|
||||
|
||||
print("A plus:")
|
||||
sp.pprint(Aplus)
|
||||
|
||||
print("A minus:")
|
||||
sp.pprint(Aminus)
|
||||
|
||||
print("A S:")
|
||||
sp.pprint(sp.simplify(Aminus - Aplus))
|
||||
|
||||
assert sp.simplify(S * D * S.inv() - A) == sp.Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
|
||||
print("B:")
|
||||
S, D = B.diagonalize()
|
||||
print("Diagonals:")
|
||||
sp.pprint(D[0, 0])
|
||||
sp.pprint(D[1, 1])
|
||||
sp.pprint(D[2, 2])
|
||||
print("S:")
|
||||
sp.pprint(S)
|
||||
print("SI:")
|
||||
sp.pprint(sp.simplify(S.inv()))
|
||||
|
||||
assert sp.simplify(S * D * S.inv() - B) == sp.Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
|
||||
m = abs(etav / eta) + abs(sp.sqrt(eta ** (3)) * sp.sqrt(g) / eta)
|
||||
|
||||
L = D + sp.Matrix([[m, 0, 0], [0, m, 0], [0, 0, m]])
|
||||
Bplus = sp.simplify(S * L * S.inv()) / 2
|
||||
L = D - sp.Matrix([[m, 0, 0], [0, m, 0], [0, 0, m]])
|
||||
Bminus = sp.simplify(S * L * S.inv()) / 2
|
||||
|
||||
print("B plus:")
|
||||
sp.pprint(Bplus)
|
||||
|
||||
print("B minus:")
|
||||
sp.pprint(Bminus)
|
||||
|
||||
print("B S:")
|
||||
sp.pprint(sp.simplify(Bminus - Bplus))
|
||||
|
||||
assert sp.simplify((Bplus + Bminus) - B) == sp.Matrix([[0, 0, 0], [0, 0, 0], [0, 0, 0]])
|
||||
|
||||
breakpoint()
|
||||
code = codegen(
|
||||
[("Aplus", Aplus), ("Aminus", Aminus), ("Bplus", Bplus), ("Bminus", Bminus)], "rust"
|
||||
)
|
||||
with open("autogen.rs", "w") as f:
|
||||
f.write(code[0][1])
|
||||
155
shallow_water/misc/report.tex
Normal file
155
shallow_water/misc/report.tex
Normal file
@@ -0,0 +1,155 @@
|
||||
\documentclass[british]{scrartcl}
|
||||
\usepackage[utf8]{inputenc}
|
||||
\usepackage[T1]{fontenc}
|
||||
|
||||
\usepackage{babel}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{amsthm}
|
||||
\usepackage{physics}
|
||||
\usepackage{cleveref}
|
||||
\usepackage{csquotes}
|
||||
|
||||
\newtheorem{definition}{Definition}
|
||||
|
||||
\title{Shallow Water Equations}
|
||||
\author{Magnus Ulimoen}
|
||||
|
||||
\newcommand{\ve}[1]{\mathbf{#1}}
|
||||
\newcommand{\kron}{\otimes}
|
||||
|
||||
|
||||
\begin{document}
|
||||
\maketitle
|
||||
|
||||
\section{Introduction}
|
||||
|
||||
This is an attempt to solve the shallow water equations (SWE) using the SBP-SAT method. Equations and assumptions are laid out in this document.
|
||||
|
||||
|
||||
\section{Formulation}
|
||||
|
||||
The SWE takes the conservative form (assuming constant density)
|
||||
\begin{gather}
|
||||
\label{eq:conservative_formulation}
|
||||
\pdv{\ve{q}}{t} + \pdv{\ve{E}}{x} + \pdv{\ve{F}}{y} = 0\\
|
||||
\ve{q} = \begin{bmatrix} \eta \\ \eta u \\ \eta v \end{bmatrix}, \hspace{1em}
|
||||
\ve{E} = \begin{bmatrix} \eta u \\ \eta u^2 + \frac{1}{2} g \eta^2 \\ \eta u v\end{bmatrix}, \hspace{1em}
|
||||
\ve{F} = \begin{bmatrix} \eta v \\ \eta u v \\ \eta v^2 + \frac{1}{2} g \eta^2 \end{bmatrix}
|
||||
\end{gather}
|
||||
where $\eta$ is height above surface, $u$ fluid motion along $x$, $v$ fluid motion along $y$, $g$ acceleration due to gravity.
|
||||
|
||||
The equation is discretised in space on a regular cartesian grid, and the operators $\pdv{x}$ is approximated using the SBP operator $D_1$ (likewise for $\pdv{y}$). This operator obeys
|
||||
\begin{definition}
|
||||
An operator $D_1$ approximating $\pdv{x}$ on the form
|
||||
\begin{gather*}
|
||||
D_1 = H^{-1} \left( Q + \frac{1}{2} B \right)
|
||||
\end{gather*}
|
||||
where
|
||||
\begin{gather*}
|
||||
Q + Q^T = 0, \\
|
||||
B = -e_0 e_0^T + e_n e_n^T, \\
|
||||
H = H^T, \\
|
||||
x^T H x > 0\,\,\, \forall x \neq 0
|
||||
\end{gather*}
|
||||
is an SBP operator
|
||||
\end{definition}
|
||||
|
||||
Applying this to a discretised version of \cref{eq:conservative_formulation}:
|
||||
\begin{gather}
|
||||
\pdv{\ve{q}}{t} + D_x \ve{E} + D_y \ve{F} = 0
|
||||
\end{gather}
|
||||
where
|
||||
\begin{gather}
|
||||
D_x = I_3 \kron I_y \kron D_1\\
|
||||
D_y = I_3 \kron D_1 \kron I_x
|
||||
\end{gather}
|
||||
and $\ve{q}$ is a \enquote{flattening} of the vector (x, y, the three fields). This formulation can be discretised in time using some appropriate scheme (eg. Runge Kutta 4).
|
||||
|
||||
\subsection{Energy stability}
|
||||
To ensure stability, we must ensure no variable grows without bounds. First the continous case, taking the inner product left and right with q (and shifting stuff):
|
||||
\begin{gather}
|
||||
\left(q, \pdv{q}{t}\right) + \left(\pdv{q}{t}, q\right)
|
||||
= - \left(q, \pdv{\ve{E}}{x} + \pdv{\ve{F}}{y}\right)
|
||||
- \left(\pdv{\ve{E}}{x} + \pdv{\ve{F}}{y}, q\right)
|
||||
= \frac{1}{2}\pdv{q^2}{t}
|
||||
\end{gather}
|
||||
Let us linearise this equation, with
|
||||
\begin{gather}
|
||||
A = \pdv{E}{q}, \hspace{1em} B = \pdv{F}{q}
|
||||
\end{gather}
|
||||
which gives us
|
||||
\begin{align}
|
||||
\nonumber \left(q, \pdv{q}{t}\right) &= - \left( q, A \pdv{q}{x} + B \pdv{q}{x} \right) \\
|
||||
\nonumber &= {\left(q_x, A q\right)} - {\left[ q^T A q \right]}_{x_0}^{x_n}
|
||||
+ {\left(q, B q_y \right)} - {\left[ q^T B q \right]}_{y_0}^{y_n} \\
|
||||
&= {\left(A^T q_x, q\right)} - {\left[ q^T A q \right]}_{x_0}^{x_n}
|
||||
+ {\left(B^T q, q_y \right)} - {\left[ q^T B q \right]}_{y_0}^{y_n}
|
||||
\end{align}
|
||||
and the following
|
||||
\begin{align}
|
||||
\nonumber \frac{1}{2}\pdv{q^2}{t} &= \left( q, \pdv{q}{t} \right) + \left( \pdv{q}{t} \right) \\
|
||||
&= (q_x, Aq) - (Aq_x, q) - \left[q^T (A + A^T) q\right] + (q_y, B q) - (B q_y, q) - [q^T(B + B^T)q]
|
||||
\end{align}
|
||||
|
||||
By symmetrising and changing variables ($\hat{q} = Sq$) we can find a suitable form which allows all the $(q, F) - (F, q)$ forms to disappear. It might not be fully symmetrisable in both $x,y$ simultaneously, but this can be skippped\ldots
|
||||
Change of coordinates can be done within the integral,
|
||||
\begin{gather}
|
||||
q^T A q_x, \hspace{1em} \hat{q} = Sq, \hspace{1em} S^T \hat{q} = q \\
|
||||
q^T A q_x = {(S^T \hat{q})}^T A (S^T \hat{q}) = \hat{q}^T S S^T D S S^T \hat{q} = \hat{q}D\hat{q}
|
||||
\end{gather}
|
||||
This does not change anything within the integral, and shows the symmetrisation necessary to make the two forms cancel. The energy is therefore bounded by the boundary terms (which are omitted in the continous form). It is enough to bound
|
||||
\begin{gather}
|
||||
A^- \text{on the right} \\
|
||||
A^+ \text{on the left} \\
|
||||
B^- \text{on the top} \\
|
||||
B^+ \text{on the bottom}
|
||||
\end{gather}
|
||||
|
||||
\subsubsection{Discrete case}
|
||||
In this section we will determine the penalty parameter $\tau$ for all directions.
|
||||
\begin{gather}
|
||||
\pdv{q}{t} = - D_x (Aq) - D_y (Aq) \\
|
||||
\end{gather}
|
||||
\begin{align}
|
||||
{\left(q, \pdv{q}{t}\right)}_H &= q^T (H_y \kron H_x) \pdv{q}{t} \\
|
||||
&= - q^T H D_x (Aq) - q^T H D_y (Aq) \\
|
||||
&= -q^T (H_y \kron H_x) (I_y \kron H^{-1} (Q + \frac{1}{2}B )) (I_y \kron A) q + q^T \\
|
||||
&- q^T (H_y \kron H_x) (H^{-1} (Q + \frac{1}{2}B ) \kron I_x) (B \kron I_x) q + q^T
|
||||
\end{align}
|
||||
|
||||
Let us just do this in one dimension, it is mostly transferable anyway\ldots
|
||||
\begin{align}
|
||||
q^T H \pdv{q}{t} &= -q^T H H^{-1} (Q + \frac{1}{2}B) A q \\
|
||||
&= - q^T (Q + \frac{1}{2} B ) A q
|
||||
\end{align}
|
||||
And the transpose
|
||||
\begin{align}
|
||||
{\left(\pdv{q}{t}\right)}^T H q = -{(Aq)}^T H D_1 q = -q^T A^T (Q + \frac{1}{2}B) q
|
||||
\end{align}
|
||||
Tranpose of this gives (can always transpose a scalar)
|
||||
\begin{align}
|
||||
q^T {(A^T (Q + \frac{1}{2}B))}^T q = q^T (Q^T + \frac{1}{2}B^T) A q
|
||||
\end{align}
|
||||
The sum of these two expressions
|
||||
\begin{gather}
|
||||
\frac{1}{2}\pdv{\norm{q}^2_H}{t} = q^T (Q + Q^T + \frac{1}{2}(B + B^T)) A q \\
|
||||
= q^T B A q = -q_0^T A q_0 + q_n^T A q_n
|
||||
\end{gather}
|
||||
We can thus control the energy rate by controlling $q_0^T A q_0$ and $q_n^T A q_n$ (that is, we set the boundaries). We do this by adding the following SAT\@.
|
||||
\begin{gather}
|
||||
SAT = \tau_0 H^{-1} e_0 e_0^T A_- (q - v) + \tau_n H^{-1} e_n e_n^T A_+ (q - v)
|
||||
\end{gather}
|
||||
Adding this to the energy form above (setting $v=0$ (data independence))
|
||||
\begin{gather}
|
||||
E = -q_0^T A q_0 + \tau_0 q_0^T A_- q_0 + \tau_0 q_0^T A_-^T q_0
|
||||
+ q_n^T A q_n + \ldots \\
|
||||
= -q_0^T (A - 2\tau_0 A_-) q_0 + \ldots
|
||||
\end{gather}
|
||||
This sets the following requirements
|
||||
\begin{gather}
|
||||
\tau_0 \ge \frac{1}{2} \\
|
||||
\tau_n \le -\frac{1}{2}
|
||||
\end{gather}
|
||||
|
||||
|
||||
\end{document}
|
||||
89
shallow_water/src/characteristics.rs
Normal file
89
shallow_water/src/characteristics.rs
Normal file
@@ -0,0 +1,89 @@
|
||||
#![allow(non_snake_case)]
|
||||
/*
|
||||
* Code generated with sympy 1.7.dev
|
||||
*
|
||||
* See http://www.sympy.org/ for more information.
|
||||
*
|
||||
* This file is part of 'project'
|
||||
*/
|
||||
use super::Float;
|
||||
|
||||
pub fn Aplus(eta: Float, etau: Float, etav: Float, g: Float) -> [[Float; 3]; 3] {
|
||||
[
|
||||
[
|
||||
(eta.powf(3.0 / 2.0) * g.sqrt() + etau.abs()) / (2.0 * eta),
|
||||
1.0 / 2.0,
|
||||
0.0,
|
||||
],
|
||||
[
|
||||
eta * g / 2.0 - etau.powi(2) / (2.0 * eta.powi(2)),
|
||||
(eta.powf(3.0 / 2.0) * g.sqrt() + 2.0 * etau + etau.abs()) / (2.0 * eta),
|
||||
0.0,
|
||||
],
|
||||
[
|
||||
-etau * etav / (2.0 * eta.powi(2)),
|
||||
etav / (2.0 * eta),
|
||||
(eta.powf(3.0 / 2.0) * g.sqrt() + etau + etau.abs()) / (2.0 * eta),
|
||||
],
|
||||
]
|
||||
}
|
||||
|
||||
pub fn Aminus(eta: Float, etau: Float, etav: Float, g: Float) -> [[Float; 3]; 3] {
|
||||
[
|
||||
[
|
||||
-(eta.powf(3.0 / 2.0) * g.sqrt() + etau.abs()) / (2.0 * eta),
|
||||
1.0 / 2.0,
|
||||
0.0,
|
||||
],
|
||||
[
|
||||
eta * g / 2.0 - etau.powi(2) / (2.0 * eta.powi(2)),
|
||||
(-eta.powf(3.0 / 2.0) * g.sqrt() + 2.0 * etau - etau.abs()) / (2.0 * eta),
|
||||
0.0,
|
||||
],
|
||||
[
|
||||
-etau * etav / (2.0 * eta.powi(2)),
|
||||
etav / (2.0 * eta),
|
||||
(-eta.powf(3.0 / 2.0) * g.sqrt() + etau - etau.abs()) / (2.0 * eta),
|
||||
],
|
||||
]
|
||||
}
|
||||
|
||||
pub fn Bplus(eta: Float, etau: Float, etav: Float, g: Float) -> [[Float; 3]; 3] {
|
||||
[
|
||||
[
|
||||
(eta.powf(3.0 / 2.0) + etav.abs()) / (2.0 * eta),
|
||||
0.0,
|
||||
1.0 / 2.0,
|
||||
],
|
||||
[
|
||||
-etau * etav / (2.0 * eta.powi(2)),
|
||||
(eta.powf(3.0 / 2.0) * g.sqrt() + etav + etav.abs()) / (2.0 * eta),
|
||||
etau / (2.0 * eta),
|
||||
],
|
||||
[
|
||||
eta * g / 2.0 - etav.powi(2) / (2.0 * eta.powi(2)),
|
||||
0.0,
|
||||
(eta.powf(3.0 / 2.0) * g.sqrt() + 2.0 * etav + etav.abs()) / (2.0 * eta),
|
||||
],
|
||||
]
|
||||
}
|
||||
|
||||
pub fn Bminus(eta: Float, etau: Float, etav: Float, g: Float) -> [[Float; 3]; 3] {
|
||||
[
|
||||
[
|
||||
-(eta.powf(3.0 / 2.0) + etav.abs()) / (2.0 * eta),
|
||||
0.0,
|
||||
1.0 / 2.0,
|
||||
],
|
||||
[
|
||||
-etau * etav / (2.0 * eta.powi(2)),
|
||||
(-eta.powf(3.0 / 2.0) * g.sqrt() + etav - etav.abs()) / (2.0 * eta),
|
||||
etau / (2.0 * eta),
|
||||
],
|
||||
[
|
||||
eta * g / 2.0 - etav.powi(2) / (2.0 * eta.powi(2)),
|
||||
0.0,
|
||||
(-eta.powf(3.0 / 2.0) * g.sqrt() + 2.0 * etav - etav.abs()) / (2.0 * eta),
|
||||
],
|
||||
]
|
||||
}
|
||||
322
shallow_water/src/lib.rs
Normal file
322
shallow_water/src/lib.rs
Normal file
@@ -0,0 +1,322 @@
|
||||
use ndarray::prelude::*;
|
||||
use sbp::*;
|
||||
|
||||
mod characteristics;
|
||||
use characteristics::{Aminus, Aplus, Bminus, Bplus};
|
||||
|
||||
const G: Float = 1.0;
|
||||
|
||||
#[derive(Clone, Debug)]
|
||||
pub struct Field(Array3<Float>);
|
||||
|
||||
impl<'a> Into<ArrayView3<'a, Float>> for &'a Field {
|
||||
fn into(self) -> ArrayView3<'a, Float> {
|
||||
self.0.view()
|
||||
}
|
||||
}
|
||||
|
||||
impl<'a> Into<ArrayViewMut3<'a, Float>> for &'a mut Field {
|
||||
fn into(self) -> ArrayViewMut3<'a, Float> {
|
||||
self.0.view_mut()
|
||||
}
|
||||
}
|
||||
|
||||
impl Field {
|
||||
fn new(ny: usize, nx: usize) -> Self {
|
||||
Self(Array3::zeros((3, ny, nx)))
|
||||
}
|
||||
fn eta(&self) -> ArrayView2<Float> {
|
||||
self.0.slice(s![0, .., ..])
|
||||
}
|
||||
fn etau(&self) -> ArrayView2<Float> {
|
||||
self.0.slice(s![1, .., ..])
|
||||
}
|
||||
fn etav(&self) -> ArrayView2<Float> {
|
||||
self.0.slice(s![2, .., ..])
|
||||
}
|
||||
/*
|
||||
fn eta_mut(&mut self) -> ArrayViewMut2<Float> {
|
||||
self.0.slice_mut(s![0, .., ..])
|
||||
}
|
||||
fn etau_mut(&mut self) -> ArrayViewMut2<Float> {
|
||||
self.0.slice_mut(s![1, .., ..])
|
||||
}
|
||||
fn etav_mut(&mut self) -> ArrayViewMut2<Float> {
|
||||
self.0.slice_mut(s![2, .., ..])
|
||||
}
|
||||
*/
|
||||
fn components_mut(
|
||||
&mut self,
|
||||
) -> (
|
||||
ArrayViewMut2<Float>,
|
||||
ArrayViewMut2<Float>,
|
||||
ArrayViewMut2<Float>,
|
||||
) {
|
||||
self.0
|
||||
.multi_slice_mut((s![0, .., ..], s![1, .., ..], s![2, .., ..]))
|
||||
}
|
||||
}
|
||||
|
||||
pub struct System {
|
||||
fnow: Field,
|
||||
fnext: Field,
|
||||
x: (Float, Float, usize),
|
||||
y: (Float, Float, usize),
|
||||
op: Box<dyn operators::UpwindOperator2d>,
|
||||
k: [Field; 4],
|
||||
}
|
||||
|
||||
impl System {
|
||||
pub fn new(x: (Float, Float, usize), y: (Float, Float, usize)) -> Self {
|
||||
let field = Field::new(y.2, x.2);
|
||||
Self {
|
||||
fnow: field.clone(),
|
||||
fnext: field.clone(),
|
||||
x,
|
||||
y,
|
||||
op: Box::new(operators::Upwind9),
|
||||
k: [field.clone(), field.clone(), field.clone(), field],
|
||||
}
|
||||
}
|
||||
pub fn nx(&self) -> usize {
|
||||
self.x.2
|
||||
}
|
||||
pub fn ny(&self) -> usize {
|
||||
self.y.2
|
||||
}
|
||||
pub fn eta(&self) -> ArrayView2<Float> {
|
||||
self.fnow.eta()
|
||||
}
|
||||
pub fn etau(&self) -> ArrayView2<Float> {
|
||||
self.fnow.etau()
|
||||
}
|
||||
pub fn etav(&self) -> ArrayView2<Float> {
|
||||
self.fnow.etav()
|
||||
}
|
||||
pub fn components_mut(
|
||||
&mut self,
|
||||
) -> (
|
||||
ArrayViewMut2<Float>,
|
||||
ArrayViewMut2<Float>,
|
||||
ArrayViewMut2<Float>,
|
||||
) {
|
||||
self.fnow.components_mut()
|
||||
}
|
||||
|
||||
fn max_dt(&self) -> Float {
|
||||
0.1 * ((self.x.1 - self.x.0) / self.x.2 as Float)
|
||||
.min((self.y.1 - self.y.0) / self.y.2 as Float)
|
||||
}
|
||||
|
||||
pub fn advance(&mut self) {
|
||||
let max_dt = self.max_dt();
|
||||
let op = &self.op;
|
||||
let rhs = move |next: &mut Field, now: &Field, _t: Float| {
|
||||
let (mut next_eta, mut next_etau, mut next_etav) = next.components_mut();
|
||||
next_eta.fill(0.0);
|
||||
next_etau.fill(0.0);
|
||||
next_etav.fill(0.0);
|
||||
|
||||
let nx = next_eta.shape()[1];
|
||||
let ny = next_eta.shape()[0];
|
||||
|
||||
if false {
|
||||
let eta = now.eta();
|
||||
for j in 0..ny {
|
||||
for i in 0..nx {
|
||||
if eta[(j, i)] <= 0.0 {
|
||||
panic!("{} {}", j, i);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
let mut temp = Array2::<Float>::zeros((ny, nx));
|
||||
|
||||
// E flux
|
||||
let mut temp_dx = temp.clone();
|
||||
azip!((dest in &mut temp, etau in now.etau()) {
|
||||
*dest = *etau;
|
||||
});
|
||||
op.diffxi(temp.view(), temp_dx.view_mut());
|
||||
next_eta.scaled_add(-1.0, &temp_dx);
|
||||
|
||||
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau()) {
|
||||
*dest = etau.powi(2)/eta + G*eta.powi(2)/2.0
|
||||
});
|
||||
op.diffxi(temp.view(), temp_dx.view_mut());
|
||||
next_etau.scaled_add(-1.0, &temp_dx);
|
||||
|
||||
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau(), etav in now.etav()) {
|
||||
*dest = etau*etav/eta
|
||||
});
|
||||
op.diffxi(temp.view(), temp_dx.view_mut());
|
||||
next_etav.scaled_add(-1.0, &temp_dx);
|
||||
|
||||
// F flux
|
||||
let mut temp_dy = temp_dx;
|
||||
azip!((dest in &mut temp, etav in now.etav()) {
|
||||
*dest = *etav;
|
||||
});
|
||||
op.diffeta(temp.view(), temp_dy.view_mut());
|
||||
next_eta.scaled_add(-1.0, &temp_dy);
|
||||
|
||||
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau(), etav in now.etav()) {
|
||||
*dest = etau*etav/eta;
|
||||
});
|
||||
op.diffeta(temp.view(), temp_dy.view_mut());
|
||||
next_etau.scaled_add(-1.0, &temp_dy);
|
||||
|
||||
azip!((dest in &mut temp, eta in now.eta(), etav in now.etav()) {
|
||||
*dest = etav.powi(2)/eta + G*eta.powi(2)/2.0
|
||||
});
|
||||
op.diffeta(temp.view(), temp_dy.view_mut());
|
||||
next_etav.scaled_add(-1.0, &temp_dy);
|
||||
|
||||
// Upwind dissipation
|
||||
if false {
|
||||
let mut temp_dx = temp_dy;
|
||||
azip!((dest in &mut temp, eta in now.eta(), etau in now.etau()) {
|
||||
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etau.abs())/eta
|
||||
});
|
||||
op.dissxi(temp.view(), temp_dx.view_mut());
|
||||
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dx) {
|
||||
*dest -= eta*diss;
|
||||
});
|
||||
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dx) {
|
||||
*dest -= etau*diss;
|
||||
});
|
||||
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dx) {
|
||||
*dest -= etav*diss;
|
||||
});
|
||||
|
||||
let mut temp_dy = temp_dx;
|
||||
azip!((dest in &mut temp, eta in now.eta(), etav in now.etav()) {
|
||||
*dest = -(eta.powf(3.0/2.0)*G.sqrt() + etav.abs())/eta
|
||||
});
|
||||
op.disseta(temp.view(), temp_dy.view_mut());
|
||||
azip!((dest in &mut next_eta, eta in now.eta(), diss in &temp_dy) {
|
||||
*dest -= eta*diss;
|
||||
});
|
||||
azip!((dest in &mut next_etau, etau in now.etau(), diss in &temp_dy) {
|
||||
*dest -= etau*diss;
|
||||
});
|
||||
azip!((dest in &mut next_etav, etav in now.etav(), diss in &temp_dy) {
|
||||
*dest -= etav*diss;
|
||||
});
|
||||
}
|
||||
|
||||
// SAT boundaries
|
||||
#[derive(Debug)]
|
||||
enum Direction {
|
||||
North,
|
||||
South,
|
||||
East,
|
||||
West,
|
||||
}
|
||||
for dir in &[
|
||||
Direction::North,
|
||||
Direction::South,
|
||||
Direction::East,
|
||||
Direction::West,
|
||||
] {
|
||||
let mut dest;
|
||||
let this;
|
||||
let other;
|
||||
match dir {
|
||||
Direction::North => {
|
||||
dest = next.0.slice_mut(s![.., ny - 1, ..]);
|
||||
this = now.0.slice(s![.., ny - 1, ..]);
|
||||
other = now.0.slice(s![.., 0_usize, ..]);
|
||||
}
|
||||
Direction::South => {
|
||||
dest = next.0.slice_mut(s![.., 0_usize, ..]);
|
||||
this = now.0.slice(s![.., 0_usize, ..]);
|
||||
other = now.0.slice(s![.., ny - 1, ..]);
|
||||
}
|
||||
Direction::East => {
|
||||
dest = next.0.slice_mut(s![.., .., nx - 1]);
|
||||
this = now.0.slice(s![.., .., nx - 1]);
|
||||
other = now.0.slice(s![.., .., 0_usize]);
|
||||
}
|
||||
Direction::West => {
|
||||
dest = next.0.slice_mut(s![.., .., 0_usize]);
|
||||
this = now.0.slice(s![.., .., 0_usize]);
|
||||
other = now.0.slice(s![.., .., nx - 1]);
|
||||
}
|
||||
}
|
||||
for ((mut dest, this), other) in dest
|
||||
.axis_iter_mut(Axis(1))
|
||||
.zip(this.axis_iter(Axis(1)))
|
||||
.zip(other.axis_iter(Axis(1)))
|
||||
{
|
||||
let tau = match dir {
|
||||
Direction::North => 1.0,
|
||||
Direction::South => -1.0,
|
||||
Direction::East => 1.0,
|
||||
Direction::West => -1.0,
|
||||
};
|
||||
let hinv = match dir {
|
||||
Direction::North | Direction::South => {
|
||||
if op.is_h2eta() {
|
||||
(ny - 2) as Float / op.heta()[0]
|
||||
} else {
|
||||
(ny - 1) as Float / op.heta()[0]
|
||||
}
|
||||
}
|
||||
Direction::East | Direction::West => {
|
||||
if op.is_h2xi() {
|
||||
(nx - 2) as Float / op.hxi()[0]
|
||||
} else {
|
||||
(nx - 1) as Float / op.hxi()[0]
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
let v = (this[0], this[1], this[2]);
|
||||
let g = (other[0], other[1], other[2]);
|
||||
let mat = match dir {
|
||||
Direction::West => Aplus(v.0, v.1, v.2, G),
|
||||
Direction::East => Aminus(v.0, v.1, v.2, G),
|
||||
Direction::North => Bminus(v.0, v.1, v.2, G),
|
||||
Direction::South => Bplus(v.0, v.1, v.2, G),
|
||||
};
|
||||
|
||||
let q = [v.0 - g.0, v.1 - g.1, v.2 - g.2];
|
||||
|
||||
let mut res = [0.0; 3];
|
||||
for j in 0..3 {
|
||||
for i in 0..3 {
|
||||
res[j] += mat[j][i] * q[i];
|
||||
}
|
||||
}
|
||||
|
||||
for i in 0..3 {
|
||||
dest[i] += tau * hinv * res[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
log::trace!("Iteration complete");
|
||||
};
|
||||
integrate::integrate::<integrate::Rk4, _, _>(
|
||||
rhs,
|
||||
&self.fnow,
|
||||
&mut self.fnext,
|
||||
&mut 0.0,
|
||||
max_dt,
|
||||
&mut self.k[..],
|
||||
);
|
||||
|
||||
std::mem::swap(&mut self.fnow, &mut self.fnext)
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_advance() {
|
||||
let mut sys = System::new((0.0, 1.0, 50), (0.0, 1.0, 50));
|
||||
sys.fnow.components_mut().0.fill(1.0);
|
||||
for _ in 0..10 {
|
||||
sys.advance();
|
||||
}
|
||||
println!("{:?}", sys.fnow);
|
||||
}
|
||||
Reference in New Issue
Block a user