added some more test including fv0, fv1 for euler

.gitignore

@@ -34,6 +34,15 @@ Makefile.in
 *.lo
 *.la
 *.tmp
+*.vtu
+*.pvtu
+__pycache__
+*.dump
+*.out
+*.log
+*.aux
+*.pickle
+*.png
 
 data
 

pydemo/CMakeLists.txt

 #dune_add_subdirs(scalar)
 #dune_add_subdirs(shallowWater)
 dune_add_subdirs(euler)
+dune_add_subdirs(chemicalreaction)

pydemo/chemicalreaction/CMakeLists.txt

+# add custom target to build tool chain for euler
+dune_python_add_test(NAME chemical_python
+                     COMMAND ${PYTHON_EXECUTABLE} test.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})

pydemo/chemicalreaction/test.py

-import time, math, sys
-from timeit import timeit
+import math
 from ufl import *
-from dune.ufl import DirichletBC, Constant
-from dune.grid import structuredGrid, cartesianDomain
-from dune.alugrid import aluCubeGrid as cubeGrid
-from dune.alugrid import aluSimplexGrid as simplexGrid
-from dune.fem import parameter
-from dune.fem.view import adaptiveLeafGridView as adaptive
-from dune.fem.space import lagrange, dgonb, raviartThomas
+from dune.ufl import DirichletBC
+from dune.grid import cartesianDomain, structuredGrid
+from dune.alugrid import aluCubeGrid
+from dune.fem.space import lagrange, dgonb
 from dune.fem.scheme import galerkin
 from dune.femdg import femDGOperator
 from dune.femdg.rk import femdgStepper
-from dune.grid import reader
 
-#gridView = cubeGrid(cartesianDomain([0,0],[1,1],[50,50]))
-domain = (reader.gmsh, "circlemeshquad.msh")
-gridView = adaptive( grid=cubeGrid(constructor=domain, dimgrid=2) )
+# gridView = structuredGrid([0,0],[2*math.pi,2*math.pi],[10,10])
+gridView = aluCubeGrid(cartesianDomain([0,0],[2*math.pi,2*math.pi],[10,10]))
 
-dimRange          = 3
-space             = dgonb( gridView, order=1, dimRange = dimRange)
-u_h               = space.interpolate(dimRange*[0], name='u_h')
-pressureSpace     = lagrange(gridView, order=1, dimRange=1)
-pressure          = pressureSpace.interpolate([1e+7],name="pressure") # Andreas: why not zero?
-velocitySpace     = raviartThomas(gridView,order=1,dimRange=1)
-transportVelocity = velocitySpace.interpolate([0,0],name="velocity")
+def computeVelocity():
+    streamSpace = lagrange(gridView, order=1, dimRange=1)
+    Psi  = streamSpace.interpolate(0,name="streamFunction")
+    u,v  = TrialFunction(streamSpace), TestFunction(streamSpace)
+    x    = SpatialCoordinate(streamSpace)
+    form = ( inner(grad(u),grad(v)) - 5*sin(x[0])*sin(x[1]) * v[0] ) * dx
+    streamScheme = galerkin([form == 0, DirichletBC(streamSpace,[0]) ])
+    streamScheme.solve(target=Psi)
+    return as_vector([-Psi[0].dx(1),Psi[0].dx(0)])
 
-# Model
 class Model:
-    dimDomain = gridView.dimension
-    dimRange = space.dimRange # three unknowns: (phi, phi cu, phi cb)^T
-
-    # Set up time
-    secperhour = 3600
-    # Andreas: what is t_start? t = Constant(t_start, "time")  # time variable
-
-    # reaction, diffusion and other coefficients
-    phi_crit = 0.1
-    rhoc = 2710  # density calcite
-    Mb = 1.0  # Molar mass suspended biomass
-    Mc = 0.1  # Molar mass calcite
-
-    Mc_rhoc = Mc / rhoc
-
-    # Well
-    Qp = 1.25e-3
-    Qw = Qp
-    # source_p = Well(mesh=mesh, Q=Qp, degree=1)  # UNCOMMENT
-    # Pressure (natural)
-    # p0_in = Constant(1.25e7)
-
-    cu_in = 2000
-    Qcu = Qp * cu_in
-
-    cb_in = 0.01
-    Qcb = Qp * cb_in
-    #
-    # Dispersivity (constant):
-    #
-    D_mol = 1.587e-9  # molecular diffusion
-    alpha_long = 0.01  # longitudinal dispersivity
-    # CONSTANT
-    # D_mech = Constant(1e-6)
-    D = D_mol
-
-    # FULL TENSOR
-    # I = Identity(dim)
-
-    kub = 0.01
-    kurease = 706.7
-    ke = kub * kurease
-    Ku = 355  # Ku = KU * rho_w = 0.355 * 1e3
-
-    # unspecific attachment rate
-    ka = 4e-6
-
-    # Biomass
-    b0 = 3.18e-7  # decay rate coeff.
-
-    # Biomass
-    b0 = 3.18e-7  # decay rate coeff.
-
-    K = 1.0e-12
-
-    # Misc. other parameters
-    mu = 0.0008  # viscosity of water
-
-    ### initial ###
-    phi0 = 0.2
-    cu0  = 0
-    cb0  = 0
-    initial = as_vector([phi0, cu0, cb0])
-
-    ### Model functions ###
-    def toPrim(U):
-        # ( phi, phi cu, phi cb ) --> (phi, cu, cb)
-        return U[0], U[1] / U[0], U[2] / U[0]
-
-    def toCons(U):
-        # ( phi, cu, cb ) --> (phi, phi cu, phi cb)
-        return U[0], U[1] * U[0], U[2] * U[0]
-
-    # circle of 0.3 around center
-    def inlet( x ):
-        return conditional(sqrt( x[0]*x[0] + x[1]*x[1] ) < 3, 1., 0. )
-
-    def darcyVelocity( p ):
-        return -1./Model.mu * Model.K * grad(p[0])
-
-    # form for pressure equation
-    def pressureForm( time ):
-        u    = TrialFunction(pressureSpace)
-        v    = TestFunction(pressureSpace)
-        x    = SpatialCoordinate(pressureSpace)
-        # n    = FacetNormal(pressureSpace)
-        dbc  = DirichletBC(pressureSpace,[ 1e7 ])
-        phi, cu, cb = Model.toPrim( u_h )
-        qu   = Model.ke * phi * Model.Mb * cb * cu / (Model.Ku + cu )
-        hour = time / Model.secperhour
-        Qw1_t = Model.inlet( x ) * conditional( hour > 20., conditional( hour < 25., Model.Qw, 0), 0 )
-        Qw2_t = Model.inlet( x ) * conditional( hour > 45., conditional( hour < 60., Model.Qw, 0), 0 )
-        pressureRhs = as_vector([ qu + Qw1_t + Qw2_t ])
-        return [ inner(grad(u),grad(v)) * dx == inner(pressureRhs, v) * dx,
-                 dbc ]
-
-    def S_i(t,x,U,DU): # or S_i for a stiff source
-        phi, cu, cb = Model.toPrim( U )
-
-        qu = Model.ke * phi * Model.Mb * cb * cu / (Model.Ku + cu )
-        qc = qu # Assumption in the report, before eq 3.12
-        qb = cb * ( phi * ( Model.b0 + Model.ka) + qc * Model.Mc_rhoc )
-
-        hour = t / Model.secperhour
-        Qu_t = Model.inlet(x) * conditional( hour > 25., conditional( hour < 45., Model.Qcu, 0), 0 )
-        Qb_t = Model.inlet(x) * conditional( hour < 20., Model.Qcb, 0 )
-        return as_vector([ -qu * Model.Mc_rhoc,
-                           -qu + Qu_t,
-                           -qb + Qb_t
-                         ])
+    transportVelocity = computeVelocity()
+    def S_e(t,x,U,DU):
+        P1 = as_vector([0.2*pi,0.2*pi]) # midpoint of first source
+        P2 = as_vector([1.8*pi,1.8*pi]) # midpoint of second source
+        f1 = conditional(dot(x-P1,x-P1) < 0.2, 1, 0)
+        f2 = conditional(dot(x-P2,x-P2) < 0.2, 1, 0)
+        f  = conditional(t<5, as_vector([f1,f2,0]), as_vector([0,0,0]))
+        r = 10*as_vector([U[0]*U[1], U[0]*U[1], -2*U[0]*U[1]])
+        return f - r
     def F_c(t,x,U):
-        phi, cu, cb = Model.toPrim(U)
-        # first flux should be zero since porosity is simply an ODE
-        return as_matrix([ Model.dimDomain*[0],
-                           [*(Model.velocity(t,x,U) * cu)],
-                           [*(Model.velocity(t,x,U) * cb)]
-                         ])
+        return as_matrix([ [*(Model.velocity(t,x,U)*u)] for u in U ])
     def maxLambda(t,x,U,n):
         return abs(dot(Model.velocity(t,x,U),n))
     def velocity(t,x,U):
-        return transportVelocity
+        return Model.transportVelocity
     def F_v(t,x,U,DU):
-        # eps * phi^(1/3) * grad U
-        return Model.D * pow(U[0], 1./3.) * DU
-    def maxDiffusion(t,x,U):
-       return Model.D * pow(U[0], 1./3.)
+        return 0.02*DU
     def physical(t,x,U):
-        #phi, cu, cb = Model.toPrim( U )
-        # U should be positive
-        #return conditional( phi>=0,1,0)*conditional(cu>=0,1,0)*conditional(cb>=0,1,0)
-        return conditional(U[0]>1e-10,1,0)*conditional(U[1]>=0,1,0)*conditional(U[2]>=0,1,0)
-    def dirichletValue(t,x,u):
-        return Model.initial
-    boundary = {1: dirichletValue}
-    endTime = 250 * secperhour
-    name = "micap"
-
-#### main program ####
-
-operator = femDGOperator(Model, space, limiter=None, threading=True)
-stepper = femdgStepper(order=1, rkType="IM") ( operator ) # Andreas: move away from 'EX'?
-# odeParam={"fem.verboserank": 0,
-#           "fem.ode.verbose": "full",
-#           "fem.solver.verbose": True,
-#           "fem.solver.newton.verbose": True,
-#           "fem.solver.method": "gmres"}
-# parameter.append( odeParam )
-
-
-for i in range(4):
-    gridView.hierarchicalGrid.globalRefine(1)
-    u_h.interpolate(Model.initial)
-    print("computing:",i,timeit("stepper(u_h,dt=360)",number=1,globals=globals()))
+        return conditional(U[0]>=0,1,0)*conditional(U[1]>=0,1,0)*\
+               conditional(U[2]>=0,1,0)
+    boundary = {range(1,5): as_vector([0,0,0])}
+
+space    = dgonb( gridView, order=3, dimRange=3)
+u_h      = space.interpolate([0,0,0], name='u_h')
+operator = femDGOperator(Model, space, limiter="scaling")
+stepper  = femdgStepper(order=3, operator=operator)
+operator.applyLimiter(u_h)
+
+# velo = expression2GF(gridView,Model.transportVelocity, order=2, name="velocity")
+vtk = gridView.sequencedVTK("chemical", subsampling=1,
+                            pointdata=[u_h],
+                            cellvector={"velocity":Model.transportVelocity})
+vtk()
+
+t        = 0
+saveStep = 0.1
+saveTime = saveStep
+while t < 0.1:
+    operator.setTime(t)
+    t += stepper(u_h)
+    if t > saveTime:
+        print(t)
+        vtk()
+        saveTime += saveStep
+vtk()

pydemo/euler/CMakeLists.txt

 # add custom target to build tool chain for euler
-dune_python_add_test(NAME euler_python
-                     COMMAND ${PYTHON_EXECUTABLE} test.py
+dune_python_add_test(NAME eulerdg_python
+                     COMMAND ${PYTHON_EXECUTABLE} testdg.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+dune_python_add_test(NAME eulerfv1_python
+                     COMMAND ${PYTHON_EXECUTABLE} testfv1.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+dune_python_add_test(NAME eulerfv0_python
+                     COMMAND ${PYTHON_EXECUTABLE} testfv0.py
                      WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})

pydemo/euler/testfv0.py

+import mpi4py.rc
+# mpi4py.rc.threaded = False
+from dune.fem import parameter
+from dune.femdg.testing import run
+
+# from euler import constant as problem
+from euler import sod as problem
+# from euler import vortex as problem
+# from euler import leVeque as problem
+# from euler import radialSod3 as problem
+
+dim = 2
+gamma = 1.4
+
+parameter.append({"fem.verboserank": -1})
+
+primitive=lambda Model,uh: {"pressure": Model.toPrim(uh)[2]}
+parameters = {"fem.ode.odesolver": "EX",
+              "fem.timeprovider.factor": 0.25,
+              "fem.ode.order": 1}
+#-----------------
+# "dgadvectionflux.method": "EULER-HLLC", "EULER-HLL", "LLF"
+# default value is 'LLF'
+#-----------------
+# femdg.limiter.tolerance: 1 (tolerance for shock indicator)
+# femdg.limiter.epsilon: 1e-8 (epsilon to avoid rounding errors)
+# femdg.limiter.admissiblefunctions:
+#    0 = only dg solution | 1 = only reconstruction | 2 = both
+#-----------------
+
+Model = problem()
+# use shorter end time for testing
+Model.endTime = 0.025
+Model.exact = None
+
+run(Model,
+    startLevel=0, limiter=None,
+    primitive=primitive, saveStep=0.16, subsamp=2,
+    dt=None,threading=False,grid="alucube", space="finitevolume",
+    parameters=parameters)
+
+# print(str(parameter))

pydemo/euler/testfv1.py

+import mpi4py.rc
+# mpi4py.rc.threaded = False
+from dune.fem import parameter
+from dune.femdg.testing import run
+
+# from euler import constant as problem
+from euler import sod as problem
+# from euler import vortex as problem
+# from euler import leVeque as problem
+# from euler import radialSod3 as problem
+
+dim = 2
+gamma = 1.4
+
+parameter.append({"fem.verboserank": -1})
+
+primitive=lambda Model,uh: {"pressure": Model.toPrim(uh)[2]}
+parameters = {"fem.ode.odesolver": "EX",
+              "fem.timeprovider.factor": 0.25,
+              "fem.ode.order": 2,
+              "femdg.limiter.admissiblefunctions": 1,
+              "femdg.limiter.tolerance": 1,
+              "femdg.limiter.epsilon": 1e-8}
+#-----------------
+# "dgadvectionflux.method": "EULER-HLLC", "EULER-HLL", "LLF"
+# default value is 'LLF'
+#-----------------
+# femdg.limiter.tolerance: 1 (tolerance for shock indicator)
+# femdg.limiter.epsilon: 1e-8 (epsilon to avoid rounding errors)
+# femdg.limiter.admissiblefunctions:
+#    0 = only dg solution | 1 = only reconstruction | 2 = both
+#-----------------
+
+Model = problem()
+# use shorter end time for testing
+Model.endTime = 0.025
+Model.exact = None
+
+run(Model,
+    startLevel=0, limiter="default",
+    primitive=primitive, saveStep=0.16, subsamp=2,
+    dt=None,threading=False,grid="alucube", space="finitevolume",
+    parameters=parameters)
+
+# print(str(parameter))

python/dune/femdg/testing.py

@@ -58,7 +58,10 @@ def run(Model, Stepper=None,
 
     # create discrete function space
     try:
-        space = create.space( space, grid, order=polOrder, dimRange=dimR)
+        if space.lower() == "finitevolume":
+            space = create.space( space, grid, dimRange=dimR)
+        else:
+            space = create.space( space, grid, order=polOrder, dimRange=dimR)
     except:
         assert space.dimRange > 0
     if modifyModel is not None: