Add source term. Add shallow water flux. Minor adjustments

README.md

@@ -7,12 +7,12 @@ Adrien Crovato, 2021
 [![License](https://img.shields.io/badge/license-Apache_2.0-blue.svg)](https://www.apache.org/licenses/LICENSE-2.0)
 
 ## Features
-dg-flo can solve one-dimensional hyperbolic partial differential equations of the form `dU/dt + dF(U)/dx = 0`, where `U` and `F(U)` are the vectors of unknowns and fluxes.  
+dg-flo can solve one-dimensional hyperbolic partial differential equations of the form `dU/dt + dF(U)/dx + S(x) = 0`, where `U`, `F(U)` and `S(x)` are the vectors of unknowns, fluxes and sources.  
 Sample problems are given under the [tests](tests/) directory for solving the
 - [x] advection equation
 - [x] Burger's equation
 - [x] Euler's equations
-- [ ] shallow water equations
+- [x] shallow water equations
 
 ## Requirements
 dg-flo needs a python 3 interpreter and its libraries, as well the `numpy` package. The `matplotlib` package is optional (needed to display the solution interactively).

num/discretization.py

@@ -39,6 +39,8 @@ def __init__(self, frm, order, flux):
         # Matrices (constants)
         self.mass = None
         self.stif = None
+        # Source term (constants)
+        self.source = None
     def __str__(self):
         return 'Discretization'
 
@@ -121,6 +123,20 @@ def __eflux(self, ifluxes, u):
             f.append(fe)
         return f
 
+    def __source(self):
+        '''Compute source term on all elements
+            since the term is constant, it is computed once and for all
+        '''
+        if not self.source:
+            self.source = []
+            if self.frm.source:
+                for e in self.elements.values():
+                    self.source.append(self.frm.source.eval(e))
+            else:
+                for e in self.elements.values():
+                    self.source.append([[0.] * len(e.evalx()) for _ in range(self.frm.nv)])
+        return self.source
+
     def compute(self, u, t):
         '''Compute rhs of equation
         '''
@@ -131,13 +147,15 @@ def compute(self, u, t):
         fi = self.__iflux(u, t)
         # Compute total fluxes on all elements
         fe = self.__eflux(fi, u)
+        # Compute sources on all elements
+        sc = self.__source()
         # Compute RHS
         rhs = np.zeros(len(u))
         i = 0
         for e in self.elements.values():
             ue = []
             for j in range(self.frm.nv):
                 ue.append(np.array(u[e.rows[j]]))
-            rhs[np.concatenate(e.rows)] = me[i].dot(-se[i].dot(np.concatenate(self.frm.flux.eval(ue))) + np.concatenate(fe[i])) # M^-1 * (-S * fp + fn)
+            rhs[np.concatenate(e.rows)] = me[i].dot(-se[i].dot(np.concatenate(self.frm.flux.eval(ue))) + np.concatenate(fe[i])) - np.concatenate(sc[i]) # M^-1 * (-S * fp + fn + M * s)
             i += 1
         return rhs

num/formulation.py

@@ -21,13 +21,14 @@
 class Formulation:
     '''Formulate a given physics
     '''
-    def __init__(self, msh, fld, nv, flux, ic, bcs):
+    def __init__(self, msh, fld, nv, flux, ic, bcs, source = None):
         # Grid and groups
         self.msh = msh # mesh
         self.field = fld # field
         # Physics
         self.nv = nv # number of variables (physical unknowns)
         self.flux = flux # flux
+        self.source = source # source term
         # Conditions
         self.ic = ic # initial condition
         self.bcs = bcs # list of boundary conditions

num/source.py

@@ -0,0 +1,37 @@
+# -*- coding: utf8 -*-
+# test encoding: à-é-è-ô-ï-€
+
+# Copyright 2021 Adrien Crovato
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+## Source term
+# Adrien Crovato
+
+class Source:
+    '''Source term
+    '''
+    def __init__(self, funs):
+        self.funs = funs # list of functions(position) for each variable
+    def __str__(self):
+        return 'Source term'
+
+    def eval(self, e):
+        '''Evaluate the source term on the nodes an element
+        '''
+        x = e.evalx()
+        s = [[0.] * len(x) for _ in range(len(self.funs))]
+        for i in range(len(x)):
+            for j in range(len(self.funs)):
+                s[j][i] = self.funs[j](x[i])
+        return s

phys/flux.py

@@ -112,7 +112,6 @@ def eval(self, u):
         f = [rho*u, rho*u^2+p, (E+p)*u]
         '''
         # Pre-pro
-        u[0] = self.__clamp(u[0], 0.01) # clamp the density
         v = u[1] / u[0] # u = rho * u / rho
         p = (self.gamma - 1) * (u[2] - 0.5 * u[1] * v) # (gamma - 1) * (E - 0.5 * rho*u*u)
         # Flux
@@ -129,7 +128,6 @@ def evald(self, u):
                   -gamma*E*u/rho + (gamma-1)*u^3, gamma*E/rho + 3*(1-gamma)/2*u^2, gamma*u]
         '''
         # Pre-pro
-        u[0] = self.__clamp(u[0], 0.01) # clamp the density
         v = u[1] / u[0] # = rho * u / rho
         e = u[2] / u[0] # = E / rho
         # Flux
@@ -143,14 +141,32 @@ def evald(self, u):
         df[2][2] = self.gamma * v
         return df
 
-    def __clamp(self, u, mn):
-        '''Clamp the solution in case in becomes non-physical
+class ShallowWater(PFlux):
+    '''Shallow water flux
+    '''
+    def __init__(self, g):
+        PFlux.__init__(self)
+        self.g = g # acceleration due to gravity
+    def __str__(self):
+        return 'Shallow water flux (g = ' + str(self.g) + ')'
+
+    def eval(self, u):
+        '''Compute the physical flux vector
+        f = [h*u, gh + u^2/2]
+        '''
+        f = [0.] * len(u)
+        f[0] = u[0] * u[1]
+        f[1] = self.g * u[0] + 0.5 * u[1] * u[1]
+        return f
+
+    def evald(self, u):
+        '''Compute the physical flux derivative matrix
+            df = [u, h;
+                  g, u]
         '''
-        # TODO EULER
-        import numpy as np
-        if isinstance(u, np.ndarray):
-            for i in range(len(u)):
-                u[i] = max([mn, u[i]])
-        else:
-            u = max([mn, u])
-        return u
+        df = [[0.] * len(u) for _ in range(len(u))]
+        df[0][0] = u[1]
+        df[0][1] = u[0]
+        df[1][0] = self.g
+        df[1][1] = u[1]
+        return df

tests/shallow.py

@@ -0,0 +1,143 @@
+# -*- coding: utf8 -*-
+# test encoding: à-é-è-ô-ï-€
+
+# Copyright 2021 Adrien Crovato
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+## Shallow water equations test
+# Adrien Crovato
+#
+# Solve the shallow water equations (solitary wave - tsunami) on a 1D grid
+
+import numpy as np
+import phys.flux as pfl
+import num.flux as nfl
+import num.source as nsrc
+import num.conditions as numc
+import num.formulation as numf
+import num.discretization as numd
+import num.tintegration as numt
+import utils.lmesh as lmsh
+import utils.writer as wrtr
+import utils.testing as tst
+
+def main(gui):
+    # Constants
+    l = 400 # domain length
+    g = 9.81 # acceleration due to gravity
+    h = [1.0, 0.01] # steady state water height
+    w = [0.1, 20] # height and width of initial wave
+    z = [20, 100] # bed slope initial and final position
+    n = 50 # number of elements
+    p = 6 # order of discretization
+    v = ['h', 'u'] # physical variables
+    cfl = 1 / (2*p+1) # half of max. Courant-Friedrichs-Levy for stability
+    # Functions
+    # bed
+    def zb(x):
+        if x < l/2 + z[0]: return 0.
+        elif x < l/2 + z[1]: return 0.5 * (h[0]-h[1]) * (1 - np.cos(np.pi/(z[1]-z[0]) * (x-l/2-z[0])))
+        else: return h[0] - h[1]
+    # source term function (bed slope * g)
+    def gdzb(x):
+        if x >= l/2 + z[0] and x < l/2 + z[1]: return g * 0.5 * np.pi * (h[0]-h[1])/(z[1]-z[0]) * np.sin(np.pi/(z[1]-z[0])*(x-l/2-z[0]))
+        else: return 0.
+    # initial and boundary conditions
+    def h0(x, t): return h[0] + w[0] * np.cos(0.5*np.pi/w[1]*(x-l/2)) - zb(x) if abs(x-l/2) < w[1] else h[0] - zb(x)
+    def u0(x, t): return 0.0
+    def hl(x, t): return h[0] - zb(x)
+    def ul(x, t): return 0.0
+    # reference solutions
+    def rfunh(x, t):
+        if x < 116:
+            return h[0]
+        elif x < 160:
+            return -0.0001*x*x + 0.0272*x - 0.7985
+        elif x < 220:
+            return h[0]
+        elif x < 272:
+            return -0.000192559021192277*x*x + 0.0810028508376352*x - 7.48640208379703
+        elif x < 300:
+            return 0.000308903323667424*x*x	- 0.186483075379383*x + 28.1505055888591
+        else:
+            return h[1]
+    def rfunu(x, t):
+        if x < 116:
+            return 0.
+        elif x < 156:
+            return 0.000353056759525554*x*x - 0.0955691190357186*x + 6.32071941182482
+        elif x < 245:
+            return -1.12067281188971e-06*x*x + 0.000289440705269957*x - 0.0167770084481066
+        elif x < 264:
+            return -3.37962520178778e-05*x*x + 0.0300418744173713*x - 5.31085043131794
+        elif x < 272:
+            return -0.00755364567412781*x*x + 4.01701626023094*x - 533.785818827200
+        else:
+            return 0.
+    funs = [rfunh, rfunu]
+    if gui:
+        gui.vars = v
+        gui.frefs = funs
+    # Parameters
+    dx = l / n # cell length
+    dt = cfl * dx / (0.5 + np.sqrt(g*h[0])) # time step
+    tmax = 20 # simulation time (l/a)
+
+    # Generate mesh
+    msh = lmsh.run(l, n)
+    fld = msh.groups[0] # field
+    inl = msh.groups[1] # inlet
+    oul = msh.groups[2] # outlet
+    # Generate formulation
+    pflx = pfl.ShallowWater(g) # physical transport flux
+    ic = numc.Initial(fld, [h0, u0]) # initial condition
+    left = numc.Boundary(inl, [numc.Dirichlet(hl), numc.Dirichlet(ul)]) # left bc
+    right = numc.Boundary(oul, [numc.Dirichlet(hl), numc.Dirichlet(ul)]) # right bc
+    formul = numf.Formulation(msh, fld, len(v), pflx, ic, [left, right], nsrc.Source([lambda x: 0., gdzb]))
+    # Generate discretization
+    nflx = nfl.LaxFried(pflx, 0.) # Lax–Friedrichs flux (0: full-upwind, 1: central)
+    disc = numd.Discretization(formul, p, nflx)
+    # Define time integration method
+    wrt = wrtr.Writer('sol', 1, v, disc)
+    tint = numt.SspRk4(disc, wrt, gui)
+    tint.run(dt, tmax)
+
+    # Test
+    uexact = [[] for _ in range(len(v))] # exact solution at element eval point
+    ucomp = [[] for _ in range(len(v))] # computed solution at element eval point
+    for c,e in disc.elements.items():
+        xe = e.evalx()
+        for j in range(len(v)):
+            for i in range(len(xe)):
+                uexact[j].append(funs[j](xe[i], tint.t))
+                ucomp[j].append(tint.u[e.rows[j][i]])
+    maxdiff = [] # infinite norm
+    nrmdiff = [] # Frobenius norm
+    for j in range(len(v)):
+        maxdiff.append(np.max(np.abs(np.array(ucomp[j]) - np.array(uexact[j]))))
+        nrmdiff.append(np.linalg.norm(np.array(ucomp[j]) - np.array(uexact[j])))
+    tests = tst.Tests()
+    tests.add(tst.Test('Max(h-h_exact)', maxdiff[0], 0., 1e-1))
+    tests.add(tst.Test('Norm(h-h_exact)', nrmdiff[0], 0., 1e-1))
+    tests.add(tst.Test('Max(u-u_exact)', maxdiff[1], 0., 1e-1))
+    tests.add(tst.Test('Norm(u-u_exact)', nrmdiff[1], 0., 1e-1))
+    tests.run()
+
+if __name__=="__main__":
+    if parse().gui:
+        import utils.gui as gui
+        main(gui.Gui())
+        input('<ENTER TO QUIT>')
+    else:
+        main(None)

utils/gui.py

@@ -83,8 +83,8 @@ def update(self, u, t, tmax):
                 lines[i*4 + 1].set_ydata(ur[v][i])
                 lines[i*4 + 2].set_ydata(us[v][i])
                 lines[i*4 + 3].set_ydata(u[e.rows[v]])
-            maxu = max([max(np.concatenate(ur[v])), max(np.concatenate(us[v]))])
-            minu = min([min(np.concatenate(ur[v])), min(np.concatenate(us[v]))])
+            maxu = max([max(np.concatenate(ur[v])), max(np.concatenate(us[v])), 0.0])
+            minu = min([min(np.concatenate(ur[v])), min(np.concatenate(us[v])), 0.0])
             if minu != maxu: ax.set_ylim(minu*1.05, maxu*1.05)
             ax.set_title('time = {:5.3f} / {:5.3f}'.format(t, tmax))
             self.figs[v].canvas.draw()