diff --git a/dune/fem-dg/models/modelwrapper.hh b/dune/fem-dg/models/modelwrapper.hh
index f28c92bbf070751e486aa53b7b7a5de26fe61437..ef2fcb9ff725e9da2161d836455b9529b604d877 100644
--- a/dune/fem-dg/models/modelwrapper.hh
+++ b/dune/fem-dg/models/modelwrapper.hh
@@ -217,6 +217,7 @@ namespace Fem
JacobianRangeType& f ) const
{
assert( hasAdvection );
+ advection_.init( local.entity() ); // TODO: this should not be required but the pulse problem fails without
advection_.diffusiveFlux( local.quadraturePoint(), u, du, f);
}
@@ -246,6 +247,7 @@ namespace Fem
RangeType& A ) const
{
assert( hasAdvection );
+ assert( 0 ); // if this is used we have to check if this is correct
// TODO: u != ubar and du != dubar
advection_.linDiffusiveFlux( u, du, local.quadraturePoint(), u, du, A);
}
@@ -292,12 +294,6 @@ namespace Fem
const bool isFluxBnd =
#endif
AdditionalType::boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
-#if 0
- std::cout << "BoundaryFlux with id=" << id << " with uLeft=" << uLeft
- << " results in flux=" << gLeft
- << " --- " << isFluxBnd
- << std::endl;
-#endif
assert( isFluxBnd );
return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
}
@@ -321,9 +317,14 @@ namespace Fem
const JacobianRangeType& jacLeft,
RangeType& gLeft ) const
{
- // TODO: need to add something to 'Addtional'?
- assert( hasDiffusion );
- return 0;
+ const DomainType normal = local.intersection().integrationOuterNormal( local.localPosition() );
+ int id = getBoundaryId( local );
+#ifndef NDEBUG
+ const bool isFluxBnd =
+#endif
+ AdditionalType::diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
+ assert( isFluxBnd );
+ return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
}
template <class LocalEvaluation>
diff --git a/dune/fem-dg/operator/dg/passtraits.hh b/dune/fem-dg/operator/dg/passtraits.hh
index f62b7dd8c551c98885531b5305a24fe6b1e8cf17..8c1ede22d01ea5dc086518d7d497eb7c9628837e 100644
--- a/dune/fem-dg/operator/dg/passtraits.hh
+++ b/dune/fem-dg/operator/dg/passtraits.hh
@@ -57,13 +57,16 @@ namespace Fem
typedef Fem::CachingQuadrature< GridPartType, 1 > FaceQuadratureType;
// Allow generalization to systems
+#if 1
typedef Fem::DiscontinuousGalerkinSpace<
FunctionSpaceType,
GridPartType, polynomialOrder,
Fem::CachingStorage > DiscreteFunctionSpaceType;
- //typedef Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType,
- // GridPartType, polOrd,
- // Fem::CachingStorage > DiscreteFunctionSpaceType;
+#else
+ typedef Fem::LegendreDiscontinuousGalerkinSpace< FunctionSpaceType,
+ GridPartType, polOrd,
+ Fem::CachingStorage > DiscreteFunctionSpaceType;
+#endif
typedef Fem::AdaptiveDiscreteFunction< DiscreteFunctionSpaceType > DestinationType;
// Indicator for Limiter
diff --git a/pydemo/euler/shock_tube.py b/pydemo/euler/shock_tube.py
index fa62fbb187795016b17606b70f47184e9f226fcd..be2122ad2f12598c7a18f4b8a9358ef9e3d148f7 100644
--- a/pydemo/euler/shock_tube.py
+++ b/pydemo/euler/shock_tube.py
@@ -73,9 +73,12 @@ def useGalerkinOp():
def useODESolver(polOrder=2, limiter='default'):
global count, t, dt, saveTime
- spaceName = "dgonb"
polOrder = polOrder
- space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
+ if False:
+ # needs change in dune/fem-dg/operator/dg/passtraits.hh
+ space = create.space("dglegendre", grid, order=polOrder, dimrange=dimR, hierarchical=False)
+ else:
+ space = create.space("dgonb", grid, order=polOrder, dimrange=dimR)
u_h = initialize(space)
rho, v, p = Model.toPrim(u_h)
operator = createFemDGSolver( Model, space, limiter=limiter )
@@ -98,7 +101,7 @@ def useODESolver(polOrder=2, limiter='default'):
t += dt
tcount += 1
if t > saveTime:
- print('dt = ', dt, 'time = ',t, 'count = ',count )
+ print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
count += 1
# rho, v, p = Model.toPrim(u_h) # is this needed - works for me # without and it should...
grid.writeVTK(name,
@@ -120,7 +123,6 @@ if True:
grid = structuredGrid(x0,x1,N)
# grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
useODESolver(2,'default') # third order with limiter
-# useODESolver(2,None) # third order with limiter
elif False:
N = [n*10 for n in N]
grid = structuredGrid(x0,x1,N)
diff --git a/pydemo/scalar/paramSolver b/pydemo/scalar/paramSolver
index 38df24aa23bb24c83587053cc7a060bc62d53c6d..08adec49888e5fe25ec9c040dbcd5bcf2125f49d 100644
--- a/pydemo/scalar/paramSolver
+++ b/pydemo/scalar/paramSolver
@@ -13,18 +13,18 @@ fem.timeprovider.factor: 0.45
fem.timeprovider.updatestep: 1
# parameter for the implicit solvers
-fem.differenceoperator.eps: 1e-12
+# fem.differenceoperator.eps: 1e-12
fem.solver.gmres.restart: 50
fem.solver.newton.verbose: true
-fem.solver.newton.linear.verbose: false
-fem.solver.verbose: false
-fem.solver.newton.maxlineariterations: 1000
-fem.solver.newton.tolerance: 1e-10
-fem.solver.newton.linabstol: 1.25e-12
-fem.solver.newton.linreduction: 1.25e-12
+fem.solver.newton.linear.verbose: true
+fem.solver.verbose: true
+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
dgdiffusionflux.method: CDG2 # diffusion flux: CDG2, CDG, BR2, IP, NIPG, BO
-dgdiffusionflux.theoryparameters: 1 # scaling with theory parameters
-dgdiffusionflux.penalty: 0.
+dgdiffusionflux.theoryparameters: 0 # scaling with theory parameters
+dgdiffusionflux.penalty: 20.
dgdiffusionflux.liftfactor: 1.0
diff --git a/pydemo/scalar/scalar.py b/pydemo/scalar/scalar.py
index 97f70c03c1f7471764713631a763021679f41be0..49690a97d87a1644847c089051679b8b08501f18 100644
--- a/pydemo/scalar/scalar.py
+++ b/pydemo/scalar/scalar.py
@@ -20,33 +20,46 @@ class Transport1D:
def jump(U,V):
return U-V
-class LinearAdvectionDiffusion1D:
- name = 'linearAD'
-
- # def F_c(t,x,U):
- # return as_matrix( [ [U[0], 0] ] )
- def F_v(t,x,U,DU):
- return 0.5*DU
- def maxLambda(t,x,U,n):
- return abs(n[0])
- def velocity(t,x,U):
- return as_vector( [0,0] )
- def physical(U):
- return 1
- def jump(U,V):
- return U-V
-def LinearAdvectionDiffusion1DDirichlet(bnd):
- class Model(LinearAdvectionDiffusion1D):
+def LinearAdvectionDiffusion1D(v,eps):
+ v = as_vector(v)
+ class Model:
+ if v is not None and eps is not None:
+ name = 'linearAD'
+ elif v is not None:
+ name = 'linearA'
+ else:
+ name = 'linearD'
+ if v is not None:
+ def F_c(t,x,U):
+ return as_matrix( [[ *(v*U[0]) ]] )
+ if eps is not None:
+ def F_v(t,x,U,DU):
+ return eps*DU
+ def maxLambda(t,x,U,n):
+ return abs(dot(v,n))
+ def velocity(t,x,U):
+ return v
+ def physical(U):
+ return 1
+ def jump(U,V):
+ return U-V
+ return Model
+def LinearAdvectionDiffusion1DMixed(v,eps,bnd):
+ class Model(LinearAdvectionDiffusion1D(v,eps)):
+ def dirichletValue(t,x,u):
+ return bnd
+ boundaryValue = {1:dirichletValue, 2:dirichletValue}
+ def zeroFlux(t,x,u,n):
+ return 0
+ diffusionBoundaryFlux = {3: zeroFlux, 4: zeroFlux}
+ return Model
+def LinearAdvectionDiffusion1DDirichlet(v,eps,bnd):
+ class Model(LinearAdvectionDiffusion1D(v,eps)):
def dirichletValue(t,x,u):
return bnd
boundaryValue = {}
for i in range(1,5): boundaryValue.update( {i: dirichletValue} )
return Model
-class LinearAdvectionDiffusion1DNeuman(LinearAdvectionDiffusion1D):
- def zeroFlux(t,x,u,n):
- return 0
- boundaryFlux = {}
- for i in range(1,5): boundaryFlux.update( {i: zeroFlux} )
class Burgers1D:
name = 'burgers'
@@ -74,13 +87,19 @@ def riemanProblem(x,x0,UL,UR):
return conditional(x<x0,as_vector(UL),as_vector(UR))
def constantTransport():
- return Transport1D, as_vector( [0.1] )
+ return Transport1D, as_vector( [0.1] ),\
+ [-1, 0], [1, 0.1], [50, 7], 0.1,\
+ "constant"
def shockTransport():
- return Transport1D, riemanProblem(x[0],-0.5, [1], [0])
+ return Transport1D, riemanProblem(x[0],-0.5, [1], [0]),\
+ [-1, 0], [1, 0.1], [50, 7], 1.0,\
+ "shockTransport"
def sinProblem():
u0 = as_vector( [sin(2*pi*x[0])] )
- return LinearAdvectionDiffusion1DDirichlet(u0), u0
+ return LinearAdvectionDiffusion1DMixed(None,0.5,u0), u0,\
+ [-1, 0], [1, 0.1], [50, 7], 0.2,\
+ "cos"
def pulse():
t = 0
@@ -89,15 +108,17 @@ def pulse():
x0 = x[0] - center[0]
x1 = x[1] - center[1]
sig2 = 0.004
- sig2PlusDt4 = sig2+(4.0*eps*t);
- xq = ( x0*cos(4.0*t) + x1*sin(4.0*t))
+ 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))
- u0 = as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
- return LinearAdvectionDiffusion1DDirichlet(u0), u0
+ ux = -4.0*x1;
+ uy = 4.0*x0;
-def cosProblem():
- return LinearAdvectionDiffusion1DNeuman, as_vector( [cos(2*pi*x[0])] )
+ u0 = as_vector( [(sig2/ (sig2PlusDt4) ) * exp (-( xq*xq + yq*yq ) / sig2PlusDt4 )] )
+ return LinearAdvectionDiffusion1DDirichlet([ux,uy],eps,u0), u0,\
+ [0, 0], [1, 1], [8,8], 1.0,\
+ "puls"
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 177b21670acc44438b2fa114b4f40c3ab2f3fc89..18673fd1dee24843496da397c628269f7fd41df7 100644
--- a/pydemo/scalar/shock_tube.py
+++ b/pydemo/scalar/shock_tube.py
@@ -7,19 +7,14 @@ from dune.femdg import createFemDGSolver
# from scalar import shockTransport as problem
# from scalar import constantTransport as probelm
-from scalar import sinProblem as problem
-# from scalar import cosProblem as problem
-# from scalar import burgersShock as problem
-# from scalar import burgersVW as problem
-# from scalar import burgersStationaryShock as problem
+# from scalar import sinProblem as problem
+from scalar import pulse as problem
-Model, initial = problem()
+Model, initial, x0,x1,N, endTime, name = problem()
parameter.append({"fem.verboserank": 0})
parameter.append("parameter")
-x0,x1,N = [-1, 0], [1, 0.1], [50, 7]
-#x0,x1,N = [0, 0], [1, 1], [8,8]
grid = structuredGrid(x0,x1,N)
# grid = create.grid("ALUSimplex", cartesianDomain(x0,x1,N))
@@ -29,7 +24,6 @@ dt = 0.001
saveStep = 0.01
saveTime = saveStep
count = 0
-endTime = 1.
def useODESolver():
global count, t, dt, saveTime
@@ -38,22 +32,22 @@ def useODESolver():
polOrder = 2
space = create.space(spaceName, grid, order=polOrder, dimrange=dimR)
u_h = space.interpolate(initial, name='u_h')
- operator = createFemDGSolver( Model, space, limiter=None )
+ operator = createFemDGSolver( Model, space, limiter=None, diffusionScheme="cdg2" )
operator.setTimeStepSize(dt)
start = time.time()
- grid.writeVTK(Model.name, celldata=[u_h], number=count)
+ grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
while t < endTime:
operator.solve(target=u_h)
dt = operator.deltaT()
t += dt
if t > saveTime:
- print('dt = ', dt, 'time = ',t, 'count = ',count )
+ print('dt = ', dt, 'time = ',t, 'count = ',count, flush=True )
count += 1
- grid.writeVTK(Model.name, celldata=[u_h], number=count)
+ grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
saveTime += saveStep
- grid.writeVTK(Model.name, celldata=[u_h], number=count)
+ grid.writeVTK(name, celldata=[u_h], number=count, subsampling=2)
print("time loop:",time.time()-start)
useODESolver()
diff --git a/python/dune/femdg/_operators.py b/python/dune/femdg/_operators.py
index b2af33d5b16130a729dd12d27acffdb3c2e8d96c..14b26b101120ed681c27478e86eff774ddc73ea7 100644
--- a/python/dune/femdg/_operators.py
+++ b/python/dune/femdg/_operators.py
@@ -149,7 +149,7 @@ def generateMethod(struct,expr, cppType, name,
# create DG operator + solver
def createFemDGSolver(Model, space,
- limiter="default"):
+ limiter="default", diffusionScheme = "cdg2"):
import dune.create as create
u = TrialFunction(space)
@@ -166,8 +166,6 @@ def createFemDGSolver(Model, space,
hasNonStiffSource = hasattr(Model,"S_ns")
if hasNonStiffSource:
advModel += inner(as_vector(S_ns(t,x,u,grad(u))),v)*dx
- else:
- advModel += inner(t*u,v)*dx
hasDiffFlux = hasattr(Model,"F_v")
if hasDiffFlux:
@@ -177,8 +175,6 @@ def createFemDGSolver(Model, space,
hasStiffSource = hasattr(Model,"S_s")
if hasStiffSource:
diffModel += inner(as_vector(S_s(t,x,u,grad(u))),v)*dx
- else:
- advModel += inner(t*u,v)*dx
advModel = create.model("elliptic",space.grid, advModel)
diffModel = create.model("elliptic",space.grid, diffModel)
@@ -270,6 +266,23 @@ 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,
+ 'bool', 'diffusionBoundaryFlux',
+ args=['const int bndId',
+ 'const double &t',
+ 'const Entity& entity', 'const Point &x',
+ 'const DomainType &normal',
+ 'const RangeType &u',
+ 'const JacobianRangeType &jac',
+ 'RangeType &result'],
+ targs=['class Entity, class Point'], static=True,
+ predefined=predefined)
boundaryValueDict = getattr(Model,"boundaryValue",None)
if boundaryValueDict is not None:
boundaryValue = {}
@@ -330,12 +343,13 @@ def createFemDGSolver(Model, space,
solverId = "Dune::Fem::Solver::Enum::fem"
formId = "Dune::Fem::Formulation::Enum::primal"
limiterId = "Dune::Fem::AdvectionLimiter::Enum::limited"
- advFluxId = "Dune::Fem::AdvectionFlux::Enum::llf"
-
+ advFluxId = "Dune::Fem::AdvectionFlux::Enum::none"
diffFluxId = "Dune::Fem::DiffusionFlux::Enum::none"
- if hasDiffFlux:
- diffFluxId = "Dune::Fem::DiffusionFlux::Enum::primal"
+ if hasDiffFlux:
+ diffFluxId = "Dune::Fem::DiffusionFlux::Enum::"+diffusionScheme
+ if hasAdvFlux:
+ advFluxId = "Dune::Fem::AdvectionFlux::Enum::llf"
if limiter == None or limiter == False or limiter.lower() == "unlimiter":
limiterId = "Dune::Fem::AdvectionLimiter::Enum::unlimited"