From c79a6a938ef7a2df8a0282d40f73de7da81d601c Mon Sep 17 00:00:00 2001
From: dedner <a.s.dedner@warwick.ac.uk>
Date: Wed, 1 Apr 2020 12:42:45 +0100
Subject: [PATCH] added some more test including fv0, fv1 for euler
---
.gitignore | 9 +
pydemo/CMakeLists.txt | 1 +
pydemo/chemicalreaction/CMakeLists.txt | 4 +
pydemo/chemicalreaction/test.py | 223 ++++++-------------------
pydemo/euler/CMakeLists.txt | 10 +-
pydemo/euler/{test.py => testdg.py} | 0
pydemo/euler/testfv0.py | 42 +++++
pydemo/euler/testfv1.py | 45 +++++
python/dune/femdg/testing.py | 5 +-
9 files changed, 168 insertions(+), 171 deletions(-)
create mode 100644 pydemo/chemicalreaction/CMakeLists.txt
rename pydemo/euler/{test.py => testdg.py} (100%)
create mode 100644 pydemo/euler/testfv0.py
create mode 100644 pydemo/euler/testfv1.py
diff --git a/.gitignore b/.gitignore
index f647b035..45b07b3e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,6 +34,15 @@ Makefile.in
*.lo
*.la
*.tmp
+*.vtu
+*.pvtu
+__pycache__
+*.dump
+*.out
+*.log
+*.aux
+*.pickle
+*.png
data
diff --git a/pydemo/CMakeLists.txt b/pydemo/CMakeLists.txt
index b9229f75..d7b63dc8 100644
--- a/pydemo/CMakeLists.txt
+++ b/pydemo/CMakeLists.txt
@@ -1,3 +1,4 @@
#dune_add_subdirs(scalar)
#dune_add_subdirs(shallowWater)
dune_add_subdirs(euler)
+dune_add_subdirs(chemicalreaction)
diff --git a/pydemo/chemicalreaction/CMakeLists.txt b/pydemo/chemicalreaction/CMakeLists.txt
new file mode 100644
index 00000000..c7605c1c
--- /dev/null
+++ b/pydemo/chemicalreaction/CMakeLists.txt
@@ -0,0 +1,4 @@
+# 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})
diff --git a/pydemo/chemicalreaction/test.py b/pydemo/chemicalreaction/test.py
index ee224fd8..4e42272b 100644
--- a/pydemo/chemicalreaction/test.py
+++ b/pydemo/chemicalreaction/test.py
@@ -1,182 +1,69 @@
-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()
diff --git a/pydemo/euler/CMakeLists.txt b/pydemo/euler/CMakeLists.txt
index f54388e7..4ae92078 100644
--- a/pydemo/euler/CMakeLists.txt
+++ b/pydemo/euler/CMakeLists.txt
@@ -1,4 +1,10 @@
# 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})
diff --git a/pydemo/euler/test.py b/pydemo/euler/testdg.py
similarity index 100%
rename from pydemo/euler/test.py
rename to pydemo/euler/testdg.py
diff --git a/pydemo/euler/testfv0.py b/pydemo/euler/testfv0.py
new file mode 100644
index 00000000..04b78d05
--- /dev/null
+++ b/pydemo/euler/testfv0.py
@@ -0,0 +1,42 @@
+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))
diff --git a/pydemo/euler/testfv1.py b/pydemo/euler/testfv1.py
new file mode 100644
index 00000000..9b186a30
--- /dev/null
+++ b/pydemo/euler/testfv1.py
@@ -0,0 +1,45 @@
+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))
diff --git a/python/dune/femdg/testing.py b/python/dune/femdg/testing.py
index 7f08c596..69f63f4b 100644
--- a/python/dune/femdg/testing.py
+++ b/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:
--
GitLab