Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • dune-fem/dune-fem-dg
  • simon.praetorius/dune-fem-dg
2 results
Show changes
Commits on Source (2)
...@@ -2,3 +2,4 @@ dune_add_subdirs(scalar) ...@@ -2,3 +2,4 @@ dune_add_subdirs(scalar)
#dune_add_subdirs(shallowWater) #dune_add_subdirs(shallowWater)
dune_add_subdirs(euler) dune_add_subdirs(euler)
dune_add_subdirs(chemicalreaction) dune_add_subdirs(chemicalreaction)
dune_add_subdirs(modeltest)
dune_python_add_test(NAME bndtest_python
SCRIPT bndtest.py
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
set_tests_properties(bndtest_python PROPERTIES TIMEOUT 2500)
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)
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)})
...@@ -12,15 +12,16 @@ from dune.fem.operator import load ...@@ -12,15 +12,16 @@ from dune.fem.operator import load
from dune.fem import parameter as parameterReader from dune.fem import parameter as parameterReader
import dune.fem 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 TypeAlias, Declaration, Variable, Struct
from dune.source.cplusplus import Method as clsMethod from dune.source.cplusplus import Method as clsMethod
from dune.source.cplusplus import SourceWriter, ListWriter, StringWriter 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 from dune.fem.utility import FemThreadPoolExecutor
# limiter can be ScalingLimiter or FV based limiter with FV type reconstructions for troubled cells # limiter can be ScalingLimiter or FV based limiter with FV type reconstructions for troubled cells
...@@ -112,8 +113,7 @@ def createOrderRedcution(domainSpace): ...@@ -112,8 +113,7 @@ def createOrderRedcution(domainSpace):
##################################################### #####################################################
## fem-dg models ## fem-dg models
##################################################### #####################################################
def femDGModels(Model, space, initialTime=0): def femDGModels(Model, space, initialTime=0, returnUFL=False):
u = TrialFunction(space) u = TrialFunction(space)
v = TestFunction(space) v = TestFunction(space)
n = FacetNormal(space) n = FacetNormal(space)
...@@ -127,7 +127,7 @@ def femDGModels(Model, space, initialTime=0): ...@@ -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 advModel = inner(t*grad(u-u),grad(v))*dx # TODO: make a better empty model
hasNonStiffSource = hasattr(Model,"S_e") hasNonStiffSource = hasattr(Model,"S_e")
if hasNonStiffSource: 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: else:
hasNonStiffSource = hasattr(Model,"S_ns") hasNonStiffSource = hasattr(Model,"S_ns")
if hasNonStiffSource: if hasNonStiffSource:
...@@ -136,19 +136,43 @@ def femDGModels(Model, space, initialTime=0): ...@@ -136,19 +136,43 @@ def femDGModels(Model, space, initialTime=0):
hasDiffFlux = hasattr(Model,"F_v") hasDiffFlux = hasattr(Model,"F_v")
if hasDiffFlux: 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: else:
diffModel = inner(t*grad(u-u),grad(v))*dx # TODO: make a better empty model diffModel = inner(t*grad(u-u),grad(v))*dx # TODO: make a better empty model
hasStiffSource = hasattr(Model,"S_i") hasStiffSource = hasattr(Model,"S_i")
if hasStiffSource: 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: else:
hasStiffSource = hasattr(Model,"S_s") hasStiffSource = hasattr(Model,"S_s")
if hasStiffSource: if hasStiffSource:
diffModel += inner(as_vector(Model.S_impl(t,x,u,grad(u))),v)*dx diffModel += inner(as_vector(Model.S_impl(t,x,u,grad(u))),v)*dx
print("Model.S_s is deprecated. Use S_i instead!") 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 conservationlaw
#from dune.fem.model import integrands as conservationlaw #from dune.fem.model import integrands as conservationlaw
# we need False here, because there are additional methods that # we need False here, because there are additional methods that
...@@ -172,6 +196,29 @@ def femDGModels(Model, space, initialTime=0): ...@@ -172,6 +196,29 @@ def femDGModels(Model, space, initialTime=0):
return [Model,advModel,diffModel] 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 ## fem-dg Operator
##################################################### #####################################################
......
from functools import reduce from functools import reduce
import ufl
from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal,\ from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal,\
Coefficient, conditional Coefficient, conditional
from dune.source.cplusplus import Variable, UnformattedExpression,\ from dune.source.cplusplus import Variable, UnformattedExpression,\
...@@ -99,10 +99,15 @@ def uflExpr(Model,space,t): ...@@ -99,10 +99,15 @@ def uflExpr(Model,space,t):
else: else:
assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk fluxes given" assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk fluxes given"
except TypeError: except TypeError:
try: if isinstance(f, ufl.core.expr.Expr):
method = f(t,x,u) # allow to pass in a ufl expression instead of a function for the boundary conditions
except TypeError: # This works for non time dependent b.c
method = f(t,x,u,n,n) 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] ) boundaryValue.update( [ (kk,method) for kk in ids] )
# add dummy values for hasBoundaryValue method # add dummy values for hasBoundaryValue method
bndReturn = True bndReturn = True
...@@ -328,10 +333,15 @@ def codeFemDg(self, *args, **kwargs): ...@@ -328,10 +333,15 @@ def codeFemDg(self, *args, **kwargs):
else: else:
assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk flux given" assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk flux given"
except TypeError: except TypeError:
try: if isinstance(f, ufl.core.expr.Expr):
method = f(t,x,u) # allow to pass in a ufl expression instead of a function for the boundary conditions
except TypeError: # This works for non time dependent b.c
method = f(t,x,u,n,n) 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] ) boundaryValue.update( [ (kk,method) for kk in ids] )
# add dummy values for hasBoundaryValue method # add dummy values for hasBoundaryValue method
bndReturn = True bndReturn = True
......