diff --git a/pydemo/euler/euler.py b/pydemo/euler/euler.py
index 7af7a857c7e7f9611e15c07e3f189f3fb8d14196..4bd50f6a5d0c04e1b7435be9ca9a0b3845425ebf 100644
--- a/pydemo/euler/euler.py
+++ b/pydemo/euler/euler.py
@@ -47,25 +47,18 @@ def CompressibleEuler(dim, gamma):
def CompressibleEulerNeuman(dim, gamma):
class Model(CompressibleEuler(dim,gamma)):
- def outflowFlux(t,x,u,n):
- return Model.F_c(t,x,u)*n
- boundaryFlux = {}
- for i in range(1,5): boundaryFlux.update( {i: outflowFlux} )
+ boundary = {range(1,5): lambda t,x,u,n: Model.F_c(t,x,u)*n}
return Model
def CompressibleEulerDirichlet(dim, gamma):
class Model(CompressibleEuler(dim,gamma)):
- def outflowValue(t,x,u):
- return u
- boundaryValue = {}
- for i in range(1,5): boundaryValue.update( {i: outflowValue} )
+ boundary = {range(1,5): lambda t,x,u: u}
return Model
def CompressibleEulerSlip(dim, gamma):
class Model(CompressibleEuler(dim,gamma)):
def outflowFlux(t,x,u,n):
_,_, p = CompressibleEuler(dim,gamma).toPrim(u)
return as_vector([ 0, *(p*n), 0 ])
- boundaryFlux = {}
- for i in range(1,5): boundaryFlux.update( {i: outflowFlux} )
+ boundary = {range(1,5): outflowFlux}
return Model
def riemanProblem(x,x0,UL,UR):
diff --git a/pydemo/euler/paramBase b/pydemo/euler/paramBase
index cba313dd29c863cffd04cb03c51fc44214167e79..5b901d39786de0382089a6b317687431116b5152 100644
--- a/pydemo/euler/paramBase
+++ b/pydemo/euler/paramBase
@@ -2,7 +2,7 @@
#------------
fem.verboserank: 0
# number of threads used in OMP program
-fem.parallel.numberofthreads: 4
+fem.parallel.numberofthreads: 1
# write diagnostics file (
# 0 don't, 1 only speedup file, 2 write all runfiles
# 3 only write 0, others at end, 4 all files at end for scaling)
diff --git a/pydemo/euler/paramSolver b/pydemo/euler/paramSolver
index 06eb1e29ace68114d1cde7094f9b9e21baa08df1..7739b84068f368d4c420717356f67beded1c5646 100644
--- a/pydemo/euler/paramSolver
+++ b/pydemo/euler/paramSolver
@@ -2,7 +2,7 @@
#---------------------
fem.ode.odesolver: EX # ode solvers: EX, IM, IMEX
-fem.ode.order: 3
+# fem.ode.order: 3
fem.ode.verbose: cfl # ode output: none, cfl, full
fem.ode.cflincrease: 1.25
fem.ode.miniterations: 95
diff --git a/pydemo/euler/shock_tube.py b/pydemo/euler/shock_tube.py
index b673a88af84c00fb590679d7fc42399a75f3d315..27ff389676d7b6a39793b70174b86aab7c23c87b 100644
--- a/pydemo/euler/shock_tube.py
+++ b/pydemo/euler/shock_tube.py
@@ -9,8 +9,9 @@ from ufl import *
gamma = 1.4
dim = 2
-from euler import sod as problem
-# from euler import radialSod3 as problem
+
+# from euler import sod as problem
+from euler import radialSod3 as problem
Model, initial, x0,x1,N, endTime, name = problem(dim,gamma)
@@ -27,6 +28,7 @@ saveStep = 0.01
saveTime = saveStep
def initialize(space):
+ return space.interpolate(initial, name='u_h')
if space.order == 0:
return space.interpolate(initial, name='u_h')
else:
@@ -83,7 +85,7 @@ def useODESolver(polOrder=2, limiter='default'):
operator = createFemDGSolver( Model, space, limiter=limiter )
# operator.setTimeStepSize(dt)
- # operator.applyLimiter( u_h );
+ operator.applyLimiter( u_h );
print("number of elements: ",grid.size(0),flush=True)
grid.writeVTK(name,
pointdata=[u_h],
@@ -101,9 +103,8 @@ def useODESolver(polOrder=2, limiter='default'):
tcount += 1
if tcount%100 == 0:
print('[',tcount,']','dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
- if False: # t > saveTime:
+ if t > saveTime:
count += 1
- # rho, v, p = Model.toPrim(u_h) # is this needed - works for me # without and it should...
grid.writeVTK(name,
pointdata=[u_h],
celldata={"pressure":p, "maxLambda":Model.maxLambda(0,0,u_h,as_vector([1,0]))},
@@ -121,6 +122,8 @@ def useODESolver(polOrder=2, limiter='default'):
if True:
# grid = structuredGrid(x0,x1,N)
grid = create.grid("ALUCube", cartesianDomain(x0,x1,N))
+ # grid.hierarchicalGrid.globalRefine(2)
+ # grid = create.view("adaptive", grid)
useODESolver(2,'default') # third order with limiter
elif False:
N = [n*10 for n in N]
diff --git a/pydemo/scalar/paramSolver b/pydemo/scalar/paramSolver
index a2df1b67f1529e73bb233ce08f54b5d8b340f950..9b8aa72fb5d7b739567744f06169ee9572942839 100644
--- a/pydemo/scalar/paramSolver
+++ b/pydemo/scalar/paramSolver
@@ -15,16 +15,17 @@ fem.timeprovider.updatestep: 1
# parameter for the implicit solvers
# fem.differenceoperator.eps: 1e-12
fem.solver.gmres.restart: 50
-fem.solver.newton.verbose: false
+fem.solver.newton.verbose: true
fem.solver.newton.linear.verbose: false
fem.solver.verbose: false
fem.solver.newton.maxlineariterations: 100000
fem.solver.newton.tolerance: 1e-07
-fem.solver.newton.linabstol: 1.25e-9
-fem.solver.newton.linreduction: 1.25e-9
+fem.solver.newton.linabstol: 2e-8
+fem.solver.newton.linreduction: 2e-8
+
dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 0 # scaling with theory parameters
-dgdiffusionflux.penalty: 20.
+dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
+dgdiffusionflux.penalty: 0.
dgdiffusionflux.liftfactor: 1.0
diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py
index 63580c72a01fc2d77ecdac8c9f5439084cd1e33f..3cb025c45fa93fc50a6101de0cae73251dc27dbb 100644
--- a/pydemo/scalar/scalar.py
+++ b/pydemo/scalar/scalar.py
@@ -6,10 +6,7 @@ class Transport1D:
def F_c(t,x,U):
return as_matrix( [ [U[0], 0] ] )
- def outflowValue(t,x,u):
- return u
- boundaryValue = {}
- for i in range(1,5): boundaryValue.update( {i: outflowValue} )
+ boundary = {range(1,5): lambda t,x,u: u}
def maxLambda(t,x,U,n):
return abs(n[0])
@@ -48,19 +45,21 @@ def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
class Model(LinearAdvectionDiffusion1D(v,eps)):
def dirichletValue(t,x,u):
return bnd(t,x)
- boundaryValue = {1:dirichletValue, 2:dirichletValue}
def zeroFlux(t,x,u,n):
return 0
- diffusionBoundaryFlux = {3: zeroFlux, 4: zeroFlux}
+ if v is not None and eps is not None:
+ boundary = {(1,2): dirichletValue, (3,4): [zeroFlux,zeroFlux] }
+ else:
+ boundary = {(1,2): dirichletValue, (3,4): zeroFlux }
return Model
def LinearAdvectionDiffusion1DDirichlet(v,eps,bnd):
class Model(LinearAdvectionDiffusion1D(v,eps)):
def dirichletValue(t,x,u):
return bnd(t,x)
- boundaryValue = {}
- for i in range(1,5): boundaryValue.update( {i: dirichletValue} )
+ boundary = { range(1,5): dirichletValue }
return Model
+# burgers problems still need to be adapted to new API
class Burgers1D:
name = 'burgers'
@@ -96,13 +95,13 @@ def shockTransport():
"shockTransport", None
def sinProblem():
- u0 = lambda t,x: as_vector( [sin(2*pi*x[0])] )
- return LinearAdvectionDiffusion1DMixed(None,0.5,u0), u0(0,x),\
+ eps = 0.5
+ u0 = lambda t,x: as_vector( [sin(2*pi*x[0])*exp(-t*eps*(2*pi)**2)] )
+ return LinearAdvectionDiffusion1DMixed(None,eps,u0), u0(0,x),\
[-1, 0], [1, 0.1], [50, 7], 0.2,\
- "cos", None
+ "cos", lambda t: u0(t,x)
-def pulse():
- eps = 0.001
+def pulse(eps=None):
center = as_vector([ 0.5,0.5 ])
x0 = x[0] - center[0]
x1 = x[1] - center[1]
@@ -112,13 +111,18 @@ def pulse():
def u0(t,x):
sig2 = 0.004
- sig2PlusDt4 = sig2+(4.0*eps*t)
+ if eps is None:
+ sig2PlusDt4 = sig2
+ else:
+ sig2PlusDt4 = sig2+(4.0*eps*t)
xq = ( x0*cos(4.0*t) + x1*sin(4.0*t)) - 0.25
yq = (-x0*sin(4.0*t) + x1*cos(4.0*t))
return as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
return LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0), u0(0,x),\
[0, 0], [1, 1], [16,16], 1.0,\
- "pulse", lambda t: u0(t,x)
+ "pulse"+str(eps), lambda t: u0(t,x)
+def diffusivePulse():
+ return pulse(0.001)
def burgersShock():
return Burgers1D, riemanProblem(x[0],-0.5,[1],[0])
diff --git a/pydemo/scalar/shock_tube.py b/pydemo/scalar/shock_tube.py
index 2ee9ad807409899d127cb7fbe5c0ac6226802ff5..9460aeafb8be10858b795ea2f5b6728b99c0c6dc 100644
--- a/pydemo/scalar/shock_tube.py
+++ b/pydemo/scalar/shock_tube.py
@@ -8,9 +8,9 @@ from dune.femdg import createFemDGSolver
from ufl import dot
# from scalar import shockTransport as problem
-# from scalar import constantTransport as probelm
# from scalar import sinProblem as problem
-from scalar import pulse as problem
+# from scalar import pulse as problem
+from scalar import diffusivePulse as problem
Model, initial, x0,x1,N, endTime, name, exact = problem()
diff --git a/python/dune/femdg/_operators.py b/python/dune/femdg/_operators.py
index 4fb348d08ce2293f2fb1b0dbb670167c911692ee..b249f6aebb9c7440e4246ba138c87de17309e3f9 100644
--- a/python/dune/femdg/_operators.py
+++ b/python/dune/femdg/_operators.py
@@ -273,13 +273,42 @@ def createFemDGSolver(Model, space,
targs=['class Intersection, class Point'], static=True,
predefined=predefined)
- boundaryFluxDict = getattr(Model,"boundaryFlux",None)
- if boundaryFluxDict is not None:
- boundaryFlux = {}
- for id,f in boundaryFluxDict.items(): boundaryFlux.update({id:f(t,x,u,n)})
- else:
- boundaryFlux = {}
- generateMethod(struct, boundaryFlux,
+
+ #####################
+ ## boundary treatment
+ boundaryDict = getattr(Model,"boundary",{})
+ boundaryAFlux = {}
+ boundaryDFlux = {}
+ boundaryValue = {}
+ for k,f in boundaryDict.items():
+ # collect all ids (could be list or range)
+ ids = []
+ try:
+ for kk in k:
+ ids += [kk]
+ except TypeError:
+ ids += [k]
+ # TODO: check that id is not used more then once
+ # figure out what type of boundary condition is used
+ if isinstance(f,tuple) or isinstance(f,list):
+ assert hasAdvFlux and hasDiffFlux, "two boundary fluxes provided for id "+str(k)+" but only one bulk flux given"
+ method = [f[0](t,x,u,n), f[1](t,x,u,n)]
+ boundaryAFlux.update( [ (kk,method[0]) for kk in ids] )
+ boundaryDFlux.update( [ (kk,method[1]) for kk in ids] )
+ else:
+ try:
+ method = f(t,x,u,n)
+ if hasAdvFlux and not hasDiffFlux:
+ boundaryAFlux.update( [ (kk,method) for kk in ids] )
+ elif not hasAdvFlux and hasDiffFlux:
+ boundaryDFlux.update( [ (kk,method) for kk in ids] )
+ else:
+ assert not (hasAdvFlux and hasDiffFlux), "one boundary fluxes provided for id "+str(k)+" but two bulk flux given"
+ except TypeError:
+ method = f(t,x,u)
+ boundaryValue.update( [ (kk,method) for kk in ids] )
+
+ generateMethod(struct, boundaryAFlux,
'bool', 'boundaryFlux',
args=['const int bndId',
'const double &t',
@@ -289,13 +318,7 @@ def createFemDGSolver(Model, space,
'RangeType &result'],
targs=['class Entity, class Point'], static=True,
predefined=predefined)
- diffusionBoundaryFluxDict = getattr(Model,"diffusionBoundaryFlux",None)
- if diffusionBoundaryFluxDict is not None:
- diffusionBoundaryFlux = {}
- for id,f in diffusionBoundaryFluxDict.items(): diffusionBoundaryFlux.update({id:f(t,x,u,n)})
- else:
- diffusionBoundaryFlux = {}
- generateMethod(struct, diffusionBoundaryFlux,
+ generateMethod(struct, boundaryDFlux,
'bool', 'diffusionBoundaryFlux',
args=['const int bndId',
'const double &t',
@@ -306,12 +329,6 @@ def createFemDGSolver(Model, space,
'RangeType &result'],
targs=['class Entity, class Point'], static=True,
predefined=predefined)
- boundaryValueDict = getattr(Model,"boundaryValue",None)
- if boundaryValueDict is not None:
- boundaryValue = {}
- for id,f in boundaryValueDict.items(): boundaryValue.update({id:f(t,x,u)})
- else:
- boundaryValue = {}
generateMethod(struct, boundaryValue,
'bool', 'boundaryValue',
args=['const int bndId',