pydemo/CMakeLists.txt

@@ -2,3 +2,4 @@ dune_add_subdirs(scalar)
 #dune_add_subdirs(shallowWater)
 dune_add_subdirs(euler)
 dune_add_subdirs(chemicalreaction)
+dune_add_subdirs(modeltest)

pydemo/modeltest/CMakeLists.txt

+dune_python_add_test(NAME bndtest_python
+                     SCRIPT bndtest.py
+                     WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
+set_tests_properties(bndtest_python PROPERTIES TIMEOUT 2500)

pydemo/modeltest/bndtest.py

+import dune.femdg
+
+import numpy as np
+from matplotlib import pyplot
+from ufl import *
+from dune.ufl import DirichletBC, Space
+from dune.grid import structuredGrid
+from dune.fem.function import gridFunction
+from dune.fem.plotting import plotPointData as plot
+from dune.fem.space import lagrange, dgonb
+from dune.fem.scheme import galerkin
+from dune.femdg import femDGUFL
+from dune.femdg.rk import femdgStepper
+
+from dune.fem import threading
+threading.use = max(4,threading.max) # use at most 4 threads
+
+gridView = structuredGrid([0,0],[1.2*np.pi,2.2*np.pi],[30,60])
+coord = SpatialCoordinate(Space(2))
+
+def test(model,Psi, show=False, tol=0):
+    if Psi is None:
+        space = lagrange(gridView, order=2, dimRange=1)
+    else:
+        space = Psi.space
+    u = space.interpolate(0,name="streamFunction")
+    # dune.generator.setNoDependencyCheck()
+    scheme = galerkin( femDGUFL(model(Psi), space) )
+    # dune.generator.setDependencyCheck()
+    scheme.solve(target=u)
+    if Psi is not None:
+        if show:
+            fig,(ax1,ax2,ax3) = pyplot.subplots(1,3,figsize=(20,5))
+            clim = None # [-1,1]
+            plot(Psi, gridView=gridView, gridLines=None,
+                      clim=clim,
+                      figure=(fig,ax1))
+            plot(u, gridView=gridView, gridLines=None,
+                    clim=clim,
+                    figure=(fig,ax2))
+            plot(Psi-u, gridView=gridView, gridLines=None, figure=(fig,ax3))
+            pyplot.show()
+
+        if tol > 0:
+            if np.max(Psi.as_numpy-u.as_numpy) >= tol:
+                print(tol,np.max(Psi.as_numpy-u.as_numpy))
+            assert np.max(Psi.as_numpy-u.as_numpy) < tol
+        if tol == 0:
+            if not np.all(np.isclose( Psi.as_numpy-u.as_numpy, 0, atol=1e-6) ):
+                print(tol,np.max(Psi.as_numpy-u.as_numpy))
+            assert np.all(np.isclose( Psi.as_numpy-u.as_numpy, 0, atol=1e-6) )
+    return u
+
+# L[u] = L_e[u] + L_i[u] = (-div(F_c) - S_e) + (-div(F_v) - S_i)
+# If       F_v = grad(u) then L[u] = -laplace(u) - S_i so S_i = -laplace(u)
+# and with F_c = bu then L[u] = -div(bu) - S_e - laplace(u) - S_i
+# Weak form: <L[u],v> = (F_c+F_v).grad(v) - (S_e+S_i)v - (F_c+F_v).n v
+class BaseModel:
+    exact = lambda x: as_vector([ sin(x[0])*cos(x[1]) ])
+    b = lambda x: as_vector([ sin(x[0])*sin(x[1]), cos(x[0])*cos(x[1]) ])
+    # Dexact   = (cos(x)cos(y), -sin(x)sin(y))
+    # lapexact = -sin(x)cos(y) - sin(x)cos(y) = -2exact
+    # b.Dexact = sin(x)sin(y) cos(x)cos(y) - cos(x)cos(y) sin(x)sin(y) #
+    #          = cos(x)sin(y) exact - cos(x)sin(y) exact
+    # top:    exact, Dexact.n = sin(x), -sin(x)sin(y) = 0  (y=2pi)
+    # bottom: exact, Dexact.n = sin(x), sin(x)sin(y)  = 0  (y=0)
+
+    def S_i(t,x,U,DU):
+        return 2 * BaseModel.exact(x)
+    def F_v(t,x,U,DU):
+        return DU
+    def S_e(t,x,U,DU):
+        return div(BaseModel.b(x))*U                # sign?
+    def F_c(t,x,U):
+        return as_vector([ BaseModel.b(x)*U[0] ])   # sign?
+
+def dirichletTest(Psi):
+    if Psi is None: # use the exact solution as Dirichlet boundary conditions
+        class Model(BaseModel):
+            boundary = {range(1,5): BaseModel.exact(coord)}
+    else: # use discrete solution as Dirichlet boundary conditions
+        class Model(BaseModel):
+            boundary = {range(1,5): Psi}
+    return Model
+
+def neumannTest(Psi,exact):
+    class Model(BaseModel):
+        def bndFlux_c(t,x,U,n):
+            if exact:
+                return dot(Model.F_c(t,x,BaseModel.exact(x)),n)
+                # - dot(Model.b(x),n) * BaseModel.exact(x)  # sign?
+            else:
+                return dot(Model.F_c(t,x,Psi),n)
+                # - dot(Model.b(x),n) * Psi                 # sign?
+        def bndFlux_v(t,x,U,DU,n):
+            if exact:
+                return dot(Model.F_v(t,x,BaseModel.exact(x),grad(BaseModel.exact(x))),n)
+                # return as_vector([ dot(grad(BaseModel.exact(x)[0]),n) ])
+            else:
+                return dot(Model.F_v(t,x,Psi,grad(Psi)),n)
+                # return as_vector([ dot(grad(Psi[0]),n) ])
+        if exact:
+            boundary = {(1,2): BaseModel.exact(coord), (3,4): (bndFlux_c,bndFlux_v)}
+        else:
+            boundary = {(1,2): Psi, (3,4): (bndFlux_c,bndFlux_v)}
+    return Model
+
+# 1) Dirichlet boundary conditions with exact solution
+Psi = test(dirichletTest,None,None)
+# 2) now use discrete solution from first step as boundary condition
+testPsi = test(dirichletTest,Psi, show=False, tol=0)
+# 3) now use Neumann conditions based on exact solution
+testPsi = test(lambda p: neumannTest(p,exact=True),Psi, show=False, tol=0)
+testPsi = test(lambda p: neumannTest(p,exact=False),Psi, show=False, tol=0.001)

pydemo/modeltest/modeltest.py

+import dune.femdg
+
+import numpy as np
+from matplotlib import pyplot
+from ufl import *
+from dune.ufl import DirichletBC, Constant
+from dune.grid import structuredGrid
+from dune.fem.function import gridFunction
+from dune.fem.plotting import plotPointData as plot
+from dune.fem.space import lagrange, dgonb
+from dune.fem.scheme import galerkin
+from dune.femdg import femDGUFL
+from dune.femdg.rk import femdgStepper
+
+from dune.fem import threading
+threading.use = max(4,threading.max) # use at most 4 threads
+
+# Given grid function u_n and
+# L[U] = L_e[U] + L_i[U] = (-div(F_c) - S_e)) + (-div(F_v) - S_i)
+# implement model for operator
+# T[U] = (U - u_n) + tau L[U]
+# So with F_v = grad(u), F_c=S_e=0:
+# T[U] = U - u_n - laplace(U) - S_i
+# so we modify S_i <= S_i - U and S_e <= S_e + u_n
+
+def timeModel(Model, old):
+    class TimeModel(Model):
+        tau = dune.ufl.Constant(0,"tau")
+        if hasattr(Model,"S_i"):
+            def S_i(t,x,U,DU): # or S_e for a non stiff source
+                return TimeModel.tau * Model.S_i(t,x,U,DU) - U
+        else:
+            def S_i(t,x,U,DU): # or S_e for a non stiff source
+                return -U
+        if hasattr(Model,"S_e"):
+            def S_e(t,x,U,DU): # or S_e for a non stiff source
+                return TimeModel.tau * Model.S_e(t,x,old,grad(old)) + old
+        else:
+            def S_e(t,x,U,DU): # or S_e for a non stiff source
+                return old
+        if hasattr(Model,"F_c"):
+            def F_c(t,x,U):
+                return TimeModel.tau * Model.F_c(t,x,old)
+        if hasattr(Model,"F_v"):
+            def F_v(t,x,U,DU):
+                return TimeModel.tau * Model.F_v(t,x,U,DU)
+    return TimeModel
+
+################################################################
+
+class VelocityModel:
+    def S_e(t,x,U,DU):
+        return conditional(x[1]<0,as_vector([3]),as_vector([-3]))
+    def F_v(t,x,U,DU):
+        return DU
+    boundary = {range(1,5): as_vector([0.])}
+
+gridView = structuredGrid([0,-np.pi],[np.pi,np.pi],[30,60])
+space = lagrange(gridView,dimRange=1)
+psi = space.interpolate(0,name="streamFunction")
+scheme = galerkin( femDGUFL(VelocityModel, space) )
+scheme.solve(target=psi)
+velocity = as_vector([-psi[0].dx(1),psi[0].dx(0)])
+# gridFunction(velocity).plot(gridLines=None, vectors=[0,1], block=False)
+# pyplot.show()
+
+#############################################################
+
+class ChemicalModel:
+    def S_e(t,x,U,DU):
+        P1 = as_vector([0.2*pi,-0.8*pi]) # midpoint of first source
+        P2 = as_vector([0.2*pi, 0.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) # issue here - other code has this as (f-r)
+    def F_c(t,x,U):
+        return as_matrix([ [*(velocity*u)] for u in U ])
+    def F_v(t,x,U,DU):
+        return 0.02*DU
+    boundary = {range(1,5): as_vector([0,0,0])}
+
+space = lagrange(gridView,dimRange=3)
+c_old = space.interpolate(as_vector([0,0,0]),name="concentrations")
+c_new = space.interpolate(as_vector([0,0,0]),name="concentrations")
+Model = timeModel(ChemicalModel, c_old)
+scheme = galerkin( femDGUFL(Model, space) )
+
+###################################
+
+tau = 0.01
+saveInterval = 100*tau
+nextSaveTime = saveInterval
+Model.tau.value = tau
+endTime = 10
+t = 0
+while t<endTime:
+    c_old.assign(c_new)
+    scheme.model.time = t
+    scheme.solve(target=c_new)
+    t += tau
+    if t > nextSaveTime:
+        print("time=",t,flush=True)
+        # c_new.plot()
+        nextSaveTime += saveInterval
+
+from dune.fem.plotting import plotComponents
+from matplotlib import ticker
+# plotComponents(c_new, gridLines=None, level=1,
+#                colorbar={"orientation":"horizontal", "ticks":ticker.MaxNLocator(nbins=4)})

python/dune/femdg/_operators.py

@@ -12,15 +12,16 @@ from dune.fem.operator import load
 from dune.fem import parameter as parameterReader
 import dune.fem
 
-from dune.ufl import Constant
+from dune.ufl import Constant, DirichletBC
 
 from dune.source.cplusplus import TypeAlias, Declaration, Variable, Struct
 from dune.source.cplusplus import Method as clsMethod
 from dune.source.cplusplus import SourceWriter, ListWriter, StringWriter
 
-from ufl import SpatialCoordinate,TestFunction,TrialFunction,as_vector,dx,grad,inner,FacetNormal
+from ufl import SpatialCoordinate,TestFunction,TrialFunction,as_vector,dx,ds, grad,inner,FacetNormal
+
+from dune.femdg.patch import transform, uflExpr
 
-from dune.femdg.patch import transform
 from dune.fem.utility import FemThreadPoolExecutor
 
 # limiter can be ScalingLimiter or FV based limiter with FV type reconstructions for troubled cells
@@ -112,8 +113,7 @@ def createOrderRedcution(domainSpace):
 #####################################################
 ## fem-dg models
 #####################################################
-def femDGModels(Model, space, initialTime=0):
-
+def femDGModels(Model, space, initialTime=0, returnUFL=False):
     u = TrialFunction(space)
     v = TestFunction(space)
     n = FacetNormal(space)
@@ -127,7 +127,7 @@ def femDGModels(Model, space, initialTime=0):
         advModel = inner(t*grad(u-u),grad(v))*dx    # TODO: make a better empty model
     hasNonStiffSource = hasattr(Model,"S_e")
     if hasNonStiffSource:
-        advModel += inner(as_vector(Model.S_e(t,x,u,grad(u))),v)*dx   # (-div F_v + S_e) * v
+        advModel += inner(as_vector(Model.S_e(t,x,u,grad(u))),v)*dx   # (-div F_c + S_e) * v
     else:
         hasNonStiffSource = hasattr(Model,"S_ns")
         if hasNonStiffSource:
@@ -136,19 +136,43 @@ def femDGModels(Model, space, initialTime=0):
 
     hasDiffFlux = hasattr(Model,"F_v")
     if hasDiffFlux:
-        diffModel = inner(Model.F_v(t,x,u,grad(u)),grad(v))*dx
+        diffModel = inner(Model.F_v(t,x,u,grad(u)),grad(v))*dx  # -div F_v v
     else:
         diffModel = inner(t*grad(u-u),grad(v))*dx   # TODO: make a better empty model
 
     hasStiffSource = hasattr(Model,"S_i")
     if hasStiffSource:
-        diffModel += inner(as_vector(Model.S_i(t,x,u,grad(u))),v)*dx
+        diffModel += inner(as_vector(Model.S_i(t,x,u,grad(u))),v)*dx # (-div F_v + S_i) v
     else:
         hasStiffSource = hasattr(Model,"S_s")
         if hasStiffSource:
             diffModel += inner(as_vector(Model.S_impl(t,x,u,grad(u))),v)*dx
             print("Model.S_s is deprecated. Use S_i instead!")
 
+    if returnUFL:
+        # need to extract boundary information from Model
+        maxWaveSpeed, velocity, diffusionTimeStep, physical, jump,\
+           boundaryAFlux, boundaryDFlux, boundaryValue, hasBoundaryValue,\
+           physicalBound = uflExpr(Model,space,t)
+        dirichletBCs = [ DirichletBC(space,item[1],item[0])
+                for item in boundaryValue.items() ]
+        boundaryDFlux = -sum( [ inner(item[1],v) * ds(item[0])
+                for item in boundaryDFlux.items() ] ) # keep all forms on left hand side
+        boundaryAFlux = -sum( [ inner(item[1],v) * ds(item[0])
+                for item in boundaryAFlux.items() ] ) # keep all forms on left hand side
+        return (advModel, diffModel,
+           {"maxWaveSpeed":maxWaveSpeed,
+            "velocity":velocity,
+            "diffusionTimeStep":diffusionTimeStep,
+            "physical":physical,
+            "dirichletBCs":dirichletBCs,
+            "boundaryAFlux":boundaryAFlux,
+            "boundaryDFlux":boundaryDFlux,
+            "boundaryValue":boundaryValue,
+            "hasBoundaryValue":hasBoundaryValue,
+            "physicalBound":physicalBound
+           })
+
     from dune.fem.model import conservationlaw
     #from dune.fem.model import integrands as conservationlaw
     # we need False here, because there are additional methods that
@@ -172,6 +196,29 @@ def femDGModels(Model, space, initialTime=0):
 
     return [Model,advModel,diffModel]
 
+def femDGUFL(Model, space, initialTime=0,*, returnFull=False):
+    class M(Model):
+        if hasattr(Model,"S_e"):
+            def S_e(t,x,U,DU):
+                return -Model.S_e(t,x,U,DU)
+        if hasattr(Model,"S_i"):
+            def S_i(t,x,U,DU):
+                return -Model.S_i(t,x,U,DU)
+        if hasattr(Model,"F_c"):
+            def F_c(t,x,U):
+                return -Model.F_c(t,x,U)
+    advModel,diffModel,otherExpr = femDGModels(M, space, initialTime, returnUFL=True)
+    otherExpr["boundaryAFlux"] = -otherExpr["boundaryAFlux"]
+    form = advModel+diffModel + otherExpr["boundaryAFlux"] + otherExpr["boundaryDFlux"]
+
+    if not returnFull:
+        return [ form == 0, *otherExpr["dirichletBCs"] ]
+    else:
+        otherExpr["advModel"] = advModel
+        otherExpr["diffModel"] = diffModel
+        otherExpr["form"] = form
+        return otherExpr
+
 #####################################################
 ## fem-dg Operator
 #####################################################

python/dune/femdg/patch.py

 
 from functools import reduce
-
+import ufl
 from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal,\
                 Coefficient, conditional
 from dune.source.cplusplus import Variable, UnformattedExpression,\
@@ -99,10 +99,15 @@ def uflExpr(Model,space,t):
                 else:
                     assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk fluxes given"
             except TypeError:
-                try:
-                    method = f(t,x,u)
-                except TypeError:
-                    method = f(t,x,u,n,n)
+                if isinstance(f, ufl.core.expr.Expr):
+                    # allow to pass in a ufl expression instead of a function for the boundary conditions
+                    # This works for non time dependent b.c
+                    method = f
+                else:
+                    try:
+                        method = f(t,x,u)
+                    except TypeError:
+                        method = f(t,x,u,n,n)
                 boundaryValue.update( [ (kk,method) for kk in ids] )
                 # add dummy values for hasBoundaryValue method
                 bndReturn = True
@@ -328,10 +333,15 @@ def codeFemDg(self, *args, **kwargs):
                 else:
                     assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk flux given"
             except TypeError:
-                try:
-                    method = f(t,x,u)
-                except TypeError:
-                    method = f(t,x,u,n,n)
+                if isinstance(f, ufl.core.expr.Expr):
+                    # allow to pass in a ufl expression instead of a function for the boundary conditions
+                    # This works for non time dependent b.c
+                    method = f
+                else: # must be suitable functions instead
+                    try:
+                        method = f(t,x,u)
+                    except TypeError:
+                        method = f(t,x,u,n,n)
                 boundaryValue.update( [ (kk,method) for kk in ids] )
                 # add dummy values for hasBoundaryValue method
                 bndReturn = True