diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0f54d49c50652d300d3e7294bab27414d0dc158d..5410c926204ebb36c37364562768e8b0b22f8a58 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -11,10 +11,10 @@ variables:
before_script:
- . /duneci/bin/duneci-init-job
- - python3 -m venv /duneci/modules/dune-pip
- - source /duneci/modules/dune-pip/bin/activate
- - pip install --upgrade pip
- - pip install ufl numpy matplotlib mpi4py portalocker
+ #- python3 -m venv /duneci/modules/dune-pip
+ #- source /duneci/modules/dune-pip/bin/activate
+ #- pip install --upgrade pip
+ #- pip install ufl numpy matplotlib mpi4py portalocker
- duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-geometry.git
- duneci-install-module https://gitlab.dune-project.org/core/dune-grid.git
@@ -24,24 +24,21 @@ before_script:
- duneci-install-module https://gitlab.dune-project.org/dune-fem/dune-fem.git
debian-11-gcc-9-17:
- image: duneci/debian:11
+ image: registry.dune-project.org/docker/ci/debian:11
script:
- - source /duneci/modules/dune-pip/bin/activate
+ #- source /duneci/modules/dune-pip/bin/activate
# issue with setup-dunepy: dune-fem-dg not yet build so dependency in dune-fem-dg fails
# - python /duneci/modules/dune-python/bin/setup-dunepy.py --opts=$CI_PROJECT_DIR/scripts/opts/ci-gcc.opts install
- duneci-standard-test
variables:
DUNECI_TOOLCHAIN: gcc-9-17
- # DUNE_LOG_FORMAT: '%(asctime)s - %(name)s - %(levelname)s - %(message)s'
- # DUNE_LOG_LEVEL: DEBUG
- # MAKEFLAGS: "--verbose --output-sync=target --print-directory"
-ubuntu:18.04--gcc:
- image: duneci/ubuntu:18.04
- script:
- - source /duneci/modules/dune-pip/bin/activate
- # issue with setup-dunepy: dune-fem-dg not yet build so dependency in dune-fem-dg fails
- # - python /duneci/modules/dune-python/bin/setup-dunepy.py --opts=$CI_PROJECT_DIR/scripts/opts/ci-gcc.opts install
- - duneci-standard-test
- variables:
- DUNECI_TOOLCHAIN: gcc-7-17
+ #ubuntu:20.04--gcc:
+ #image: duneci/ubuntu:18.04
+ #script:
+ #- source /duneci/modules/dune-pip/bin/activate
+ # # issue with setup-dunepy: dune-fem-dg not yet build so dependency in dune-fem-dg fails
+ # # - python /duneci/modules/dune-python/bin/setup-dunepy.py --opts=$CI_PROJECT_DIR/scripts/opts/ci-gcc.opts install
+ #- duneci-standard-test
+ #variables:
+ #DUNECI_TOOLCHAIN: gcc-7-17
diff --git a/dune/fem-dg/algorithm/sub/evolution.hh b/dune/fem-dg/algorithm/sub/evolution.hh
index 6ca3864cc9602a572528d5270956bc88b8eba6f8..054949055a5456fc5a0c45021d33d421eca1d4aa 100644
--- a/dune/fem-dg/algorithm/sub/evolution.hh
+++ b/dune/fem-dg/algorithm/sub/evolution.hh
@@ -296,6 +296,20 @@ namespace Fem
diagnostics_.registerData( "OverallTimer", &overallTime_ );
diagnostics_.registerData( "NumberOfElements", &odeSolverMonitor_.numberOfElements_ );
}
+
+#ifdef BASEFUNCTIONSET_CODEGEN_GENERATE
+ {
+ const int polOrder = solution().space().order();
+ const int maxOrder = 2 * polOrder + 1 ;
+ const int dim = GridType :: dimension;
+ std::vector< int > quadOrder( maxOrder );
+ for( int i=2; i<maxOrder+2; ++i )
+ quadOrder.push_back( i );
+ std::string filename( "autogeneratedcode_"+std::to_string( dim ) + "d_k" + std::to_string( polOrder ) + ".hh" );
+ Dune::Fem::generateCode( solution().space(), quadOrder, quadOrder, "./", filename );
+
+ }
+#endif
}
virtual void doSolve ( const int loop, TimeProviderType& tp ) override
diff --git a/dune/fem-dg/examples/advdiff/algorithmcreator.hh b/dune/fem-dg/examples/advdiff/algorithmcreator.hh
index e5d4ca8e85ddf99138afa50d00fb2199b0c522a3..0aab62457f547439ee6ff4fb2e6c2ff29690aa77 100644
--- a/dune/fem-dg/examples/advdiff/algorithmcreator.hh
+++ b/dune/fem-dg/examples/advdiff/algorithmcreator.hh
@@ -56,7 +56,10 @@ namespace Fem
typedef AlgorithmConfigurator< GridSelectorGridType,
Galerkin::Enum::default_,
Adaptivity::Enum::default_,
- DiscreteFunctionSpaces::Enum::default_, //legendre
+ //DiscreteFunctionSpaces::Enum::gausslobatto, //legendre
+ //DiscreteFunctionSpaces::Enum::gausslegendre, //legendre
+ DiscreteFunctionSpaces::Enum::legendre,
+ //DiscreteFunctionSpaces::Enum::default_,
Solver::Enum::default_,
AdvectionLimiter::Enum::default_,
Matrix::Enum::default_,
diff --git a/dune/fem-dg/examples/advdiff/models.hh b/dune/fem-dg/examples/advdiff/models.hh
index d69f9a8468da58f3b66b8b2a6cc2c830a7c7cf10..f9f2df0cf0dc9736f3a57313483fb08e9c73f884 100644
--- a/dune/fem-dg/examples/advdiff/models.hh
+++ b/dune/fem-dg/examples/advdiff/models.hh
@@ -103,8 +103,6 @@ namespace Fem
using BaseType::time_;
using BaseType::time;
- HeatEqnModel(const HeatEqnModel& otehr);
- const HeatEqnModel &operator=(const HeatEqnModel &other);
public:
/**
* \brief Constructor
diff --git a/dune/fem-dg/examples/advdiff/test/CMakeLists.txt b/dune/fem-dg/examples/advdiff/test/CMakeLists.txt
index 93fbf1c8a06a968b57b1ece08518a63228dba330..7fb9a636f2e779880b211d2929acb616ede55e6e 100644
--- a/dune/fem-dg/examples/advdiff/test/CMakeLists.txt
+++ b/dune/fem-dg/examples/advdiff/test/CMakeLists.txt
@@ -1,13 +1,13 @@
include_directories("${CMAKE_SOURCE_DIR}/dune/fem-dg/examples/advdiff/")
# specify type of grid implementation, dimension and polynomial order
-set( GRIDTYPE ALUGRID_CUBE )
-#set( GRIDTYPE SPGRID )
+#set( GRIDTYPE ALUGRID_CUBE )
+set( GRIDTYPE SPGRID )
set( GRIDDIM 2 )
-set( POLORDER 2 )
+set( POLORDER 4 )
add_definitions( "-DGRIDDIM=${GRIDDIM}" )
-add_definitions( "-DDIMRANGE=5" )
+add_definitions( "-DDIMRANGE=3" )
dune_add_test_case( NAME advdiff
SOURCES ../main.cc
diff --git a/dune/fem-dg/examples/advdiff/test/parameters/advdiff b/dune/fem-dg/examples/advdiff/test/parameters/advdiff
index fb3a7427fe8e27ec957be95fec4df826bb987f8c..a6390f72c02bb4696a463ded666b390b2cc011ba 100644
--- a/dune/fem-dg/examples/advdiff/test/parameters/advdiff
+++ b/dune/fem-dg/examples/advdiff/test/parameters/advdiff
@@ -15,7 +15,15 @@ gridsol.filename: heat-checkpoint
# GENERAL
#--------
-# femdg.stepper.printcount: 10
+fem.adaptation.method: none # none | generic | callback
+
+# coarsest level that should be present
+fem.adaptation.coarsestLevel: 3
+# finest level that should be present
+fem.adaptation.finestLevel: 6
+
+fem.eoc.steps: 1
+femdg.stepper.printcount: 100
paramfile: ../../parameter/paramBase
# PROBLEM SETUP
@@ -24,7 +32,6 @@ paramfile: ../../parameter/paramBase
# problem: heat, quasi, plaplace
problem: pulse
-fem.eoc.steps: 3
femdg.stepper.endtime: 0.15
femdg.stepper.diffusiontimestep: 1
epsilon: 0.001
diff --git a/dune/fem-dg/examples/euler/algorithmcreator.hh b/dune/fem-dg/examples/euler/algorithmcreator.hh
index 9cbb4a9674e61dadb81e432972cde4344873769a..a2822b912793893e36846f9ac4a4525d6643e1b3 100644
--- a/dune/fem-dg/examples/euler/algorithmcreator.hh
+++ b/dune/fem-dg/examples/euler/algorithmcreator.hh
@@ -62,7 +62,11 @@ namespace Fem
typedef AlgorithmConfigurator< GridSelectorGridType,
Galerkin::Enum::dg,
Adaptivity::Enum::yes,
- DiscreteFunctionSpaces::Enum::orthonormal,
+ //DiscreteFunctionSpaces::Enum::orthonormal,
+ //DiscreteFunctionSpaces::Enum::lagrange,
+ DiscreteFunctionSpaces::Enum::gausslobatto,
+ //DiscreteFunctionSpaces::Enum::gausslegendre,
+ //DiscreteFunctionSpaces::Enum::hierarchic_legendre,
Solver::Enum::fem,
AdvectionLimiter::Enum::limited,
Matrix::Enum::matrixfree,
diff --git a/dune/fem-dg/examples/euler/test/CMakeLists.txt b/dune/fem-dg/examples/euler/test/CMakeLists.txt
index 73005fd61a0aeca37523461dd51e16d027faefa3..61c78881f60a7ce2f28be904f2a25fd2d2797a0f 100644
--- a/dune/fem-dg/examples/euler/test/CMakeLists.txt
+++ b/dune/fem-dg/examples/euler/test/CMakeLists.txt
@@ -5,7 +5,7 @@ set( GRIDTYPE ALUGRID_CUBE )
#set( GRIDTYPE SPGRID_COUNT_FLOPS )
#set( GRIDTYPE SPGRID )
set( GRIDDIM 2 )
-set( POLORDER 2 )
+set( POLORDER 4 )
add_definitions( "-D${GRIDTYPE}" )
add_definitions( "-DGRIDDIM=${GRIDDIM}" )
diff --git a/dune/fem-dg/examples/euler/test/parameters/euler b/dune/fem-dg/examples/euler/test/parameters/euler
index c9d661f379929602f730fc6d68b0c33639ddcccc..e0ed5802faaecb31175c5c8fb5e690964f6cace8 100644
--- a/dune/fem-dg/examples/euler/test/parameters/euler
+++ b/dune/fem-dg/examples/euler/test/parameters/euler
@@ -8,7 +8,7 @@ problem: sod
# LIMITER SETTINGS
#-----------------
# 0 = only dg solution | 1 = only reconstruction | 2 = both
-femdg.limiter.admissiblefunctions: 1
+#femdg.limiter.admissiblefunctions: 1
# tolerance for shock indicator
femdg.limiter.tolerance: 1
# threshold for avoiding over-excessive limitation
@@ -19,7 +19,7 @@ femdg.limiter.indicatoroutput: true
# DATA WRITER
#------------
fem.io.datafileprefix: RC
-fem.io.savestep: 0.16 # SaveStep (write data every `saveStep' time period, <=0 deactivates)
+fem.io.savestep: 0.01 # SaveStep (write data every `saveStep' time period, <=0 deactivates)
fem.io.savecount: -1 # SaveCount (write data every saveCount time steps, <=0 deactivates)
@@ -29,8 +29,11 @@ gridsol.savestep: 0.5
gridsol.firstwrite: 0.5
gridsol.filename: straka-checkpoint
-fem.eoc.steps: 3
+fem.eoc.steps: 4
+femdg.stepper.printcount: 100
+# coarsest level that should be present (also initial refinement)
+fem.adaptation.coarsestLevel: 0
# STEPPER
#--------
@@ -65,13 +68,14 @@ r: 0.25 # radius of perturbation ball
# coarsest level that should be present
fem.adaptation.coarsestLevel: 0
# finest level that should be present
-fem.adaptation.finestLevel: 0
+fem.adaptation.finestLevel: 6
# DOMAIN SETUP
#-------------
fem.io.macroGridFile_1d: ../../grids/unitcube1.dgf
-fem.io.macroGridFile_2d: ../../grids/grid2d_str1d.dgf
+#fem.io.macroGridFile_2d: ../../grids/grid2d_str1d.dgf
+fem.io.macroGridFile_2d: ../../grids/grid2d_nonaffine.dgf
fem.io.macroGridFile_3d: ../../grids/unitcube3.dgf
diff --git a/dune/fem-dg/examples/euler/test/parameters/parameter b/dune/fem-dg/examples/euler/test/parameters/parameter
index b5019cedc020a1cd71becbcc844d0bcf2c7771b2..d453d8d0ab311a0919c52940130f7f965a444fca 100644
--- a/dune/fem-dg/examples/euler/test/parameters/parameter
+++ b/dune/fem-dg/examples/euler/test/parameters/parameter
@@ -1,19 +1,19 @@
-# PROBLEM SELECTION
+# PROBLEM SELECTION
#------------------
# possible: "sod" , "withman", "withmansmooth", "smooth1d" , "ffs" , "diffraction" , "shockbubble"
-problem: sod
+problem: sod
# problemflag: 0
# riemanndata: [1., -1., 1., 1., 1., 1.]
# LIMITER SETTINGS
#-----------------
-# 0 = only dg solution | 1 = only reconstruction | 2 = both
-femdg.limiter.admissiblefunctions: 1
-# tolerance for shock indicator
+# 0 = only dg solution | 1 = only reconstruction | 2 = both
+femdg.limiter.admissiblefunctions: 1
+# tolerance for shock indicator
femdg.limiter.tolerance: 1
-# threshold for avoiding over-excessive limitation
+# threshold for avoiding over-excessive limitation
femdg.limiter.limiteps: 1e-8
-# add indicator to outputvariables
+# add indicator to outputvariables
femdg.limiter.indicatoroutput: true
# DATA WRITER
@@ -23,7 +23,7 @@ fem.io.savestep: 0.01 # SaveStep (write data every `saveStep' time period, <=0 d
fem.io.savecount: -1 # SaveCount (write data every saveCount time steps, <=0 deactivates)
-# GRID SOLUTION
+# GRID SOLUTION
#--------------
gridsol.savestep: 0.5
gridsol.firstwrite: 0.5
@@ -32,12 +32,17 @@ gridsol.filename: straka-checkpoint
# GENERAL
#--------
+# coarsest level that should be present
+fem.adaptation.coarsestLevel: 3
+# finest level that should be present
+fem.adaptation.finestLevel: 3
+
paramfile: ../../parameter/paramBase
# STEPPER
#--------
-fem.eoc.steps: 4
+fem.eoc.steps: 1
femdg.stepper.starttime: 0.
femdg.stepper.endtime: 0.15
femdg.stepper.maxtimestep: 0.1
diff --git a/dune/fem-dg/examples/grids/grid2d_nonaffine.dgf b/dune/fem-dg/examples/grids/grid2d_nonaffine.dgf
new file mode 100644
index 0000000000000000000000000000000000000000..b6a197c7c60a2f971bbba38b5d8fcad878dce985
--- /dev/null
+++ b/dune/fem-dg/examples/grids/grid2d_nonaffine.dgf
@@ -0,0 +1,76 @@
+DGF
+% Elements = 16 | Vertices = 27
+
+VERTEX
+1.2500000000000000e-01 0.0000000000000000e+00
+2.5000000000000000e-01 0.0000000000000000e+00
+1.2500000000000000e-01 1.5000000000000000e-01
+2.8000000000000000e-01 1.1000000000000000e-01
+0.0000000000000000e+00 0.0000000000000000e+00
+0.0000000000000000e+00 1.2500000000000000e-01
+3.7500000000000000e-01 0.0000000000000000e+00
+5.0000000000000000e-01 0.0000000000000000e+00
+3.4500000000000000e-01 1.4500000000000000e-01
+5.0000000000000000e-01 1.2500000000000000e-01
+6.2500000000000000e-01 0.0000000000000000e+00
+5.9000000000000000e-01 1.2000000000000000e-01
+7.5000000000000000e-01 0.0000000000000000e+00
+7.2000000000000000e-01 1.5000000000000000e-01
+8.7500000000000000e-01 0.0000000000000000e+00
+1.0000000000000000e+00 0.0000000000000000e+00
+8.7500000000000000e-01 1.1500000000000000e-01
+1.0000000000000000e+00 1.2500000000000000e-01
+5.0000000000000000e-01 2.5000000000000000e-01
+6.2500000000000000e-01 2.5000000000000000e-01
+7.5000000000000000e-01 2.5000000000000000e-01
+8.7500000000000000e-01 2.5000000000000000e-01
+1.0000000000000000e+00 2.5000000000000000e-01
+3.7500000000000000e-01 2.5000000000000000e-01
+2.5000000000000000e-01 2.5000000000000000e-01
+1.2500000000000000e-01 2.5000000000000000e-01
+0.0000000000000000e+00 2.5000000000000000e-01
+#
+
+CUBE
+0 1 2 3
+4 0 5 2
+6 7 8 9
+1 6 3 8
+7 10 9 11
+10 12 11 13
+14 15 16 17
+12 14 13 16
+9 11 18 19
+11 13 19 20
+16 17 21 22
+13 16 20 21
+8 9 23 18
+3 8 24 23
+2 3 25 24
+5 2 26 25
+#
+
+BOUNDARYSEGMENTS
+2 0 1
+2 4 5
+2 4 0
+2 6 7
+2 1 6
+2 7 10
+2 10 12
+2 15 17
+2 14 15
+2 12 14
+2 18 19
+2 19 20
+2 17 22
+2 21 22
+2 20 21
+2 23 18
+2 24 23
+2 25 24
+2 5 26
+2 26 25
+#
+
+#
diff --git a/dune/fem-dg/examples/parameter/paramBase b/dune/fem-dg/examples/parameter/paramBase
index d4ae230c80c9de82a34329b3942ff8ad078d0ce1..85894f3937b58d793be647f453b2af73213d7f3c 100644
--- a/dune/fem-dg/examples/parameter/paramBase
+++ b/dune/fem-dg/examples/parameter/paramBase
@@ -5,7 +5,7 @@ fem.verboserank: 0
# OMP THREADS
#------------
# 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)
@@ -39,7 +39,7 @@ gridsol.usegridsolution: false
# STEPPER
#-----------------
-fem.eoc.steps: 1
+fem.eoc.steps: 4
femdg.stepper.printcount: 10
# maximal number of timesteps done (for debugging)
# femdg.stepper.maximaltimesteps: 10
diff --git a/dune/fem-dg/examples/parameter/paramSolver b/dune/fem-dg/examples/parameter/paramSolver
index 9c63a9fca182da81a2af9ec975fad29b59137f9a..a1623f22afb21125509f79434f3043cfdbe17452 100644
--- a/dune/fem-dg/examples/parameter/paramSolver
+++ b/dune/fem-dg/examples/parameter/paramSolver
@@ -9,7 +9,7 @@ fem.ode.miniterations: 95
fem.ode.maxiterations: 105
fem.ode.cflStart: 1.
#fem.ode.cflMax: 5
-fem.timeprovider.factor: 0.45
+fem.timeprovider.factor: 0.15
fem.timeprovider.updatestep: 1
# parameter for the implicit solvers
diff --git a/dune/fem-dg/examples/poisson/test/parameters/parameter b/dune/fem-dg/examples/poisson/test/parameters/parameter
index bb177d8eb7529097a5be48b02d7836498b5af79c..bdb9a0ab018a43ede088201926a4903c8d37fce3 100644
--- a/dune/fem-dg/examples/poisson/test/parameters/parameter
+++ b/dune/fem-dg/examples/poisson/test/parameters/parameter
@@ -1,4 +1,6 @@
-# prefix data data files
+fem.io.verboserank: 0
+
+# prefix data data files
fem.io.datafileprefix: solution
# if you want to add time stemp to eoc file name
# # to avoid overwriting the eoc files
@@ -19,7 +21,7 @@ zvelocity: 0. # the only advection part for the linear heat eqn
fem.solver.verbose: 1
fem.solver.newton.verbose: 1
fem.ode.verbose: full
-
+
# macro grid file
fem.io.macroGridFile_2d: ../../grids/unitcube2.dgf
#fem.io.macroGridFile_2d:../../grids/square2d.dgf
@@ -31,11 +33,11 @@ fem.io.macroGridFile_2d: ../../grids/unitcube2.dgf
nonConformOrigin: false
# choises are: CDG2, CDG, IP, NIPG, BO, BR2
-dgdiffusionflux.method: CDG2
+dgdiffusionflux.method: CDG2
dgdiffusionflux.theoryparameters: 1
dgdiffusionflux.lifting: id_id
#dgdiffusionflux.upwind: -1 -0.001
-#dgdiffusionflux.upwind: 0 0
+#dgdiffusionflux.upwind: 0 0
dgdiffusionflux.penalty: 0 # for CDG, CDG2, LDG
dgdiffusionflux.liftfactor: 0
diff --git a/dune/fem-dg/misc/algorithmcreatorenums.hh b/dune/fem-dg/misc/algorithmcreatorenums.hh
index e131cc4f0f93291333fb514594fb4e476c1dfc8b..bae2d3f72724c516f44c6a94ec78ef0bd2416502 100644
--- a/dune/fem-dg/misc/algorithmcreatorenums.hh
+++ b/dune/fem-dg/misc/algorithmcreatorenums.hh
@@ -77,7 +77,11 @@ namespace Fem
//! Discrete function space with Lagrange Finite Elements
lagrange,
//! p-adaptive space from dune-fem, implementing dg and lagrange
- padaptive
+ padaptive,
+ //! Lagrange space with GaussLobatto interpolation points
+ gausslobatto,
+ //! Lagrange space with GaussLegendre interpolation points
+ gausslegendre
};
}
diff --git a/dune/fem-dg/misc/algorithmcreatorselector.hh b/dune/fem-dg/misc/algorithmcreatorselector.hh
index 7939eb4f5d8782ede372e047e9f9fdc40c4e6caf..d758b46154de6a268c3cc1dc006b77ab3035ac77 100644
--- a/dune/fem-dg/misc/algorithmcreatorselector.hh
+++ b/dune/fem-dg/misc/algorithmcreatorselector.hh
@@ -51,6 +51,7 @@
#include <dune/fem-dg/assemble/primalmatrix.hh>
#include <dune/fem/space/discontinuousgalerkin.hh>
+#include <dune/fem/space/lagrange.hh>
#include <dune/fem-dg/operator/dg/operatortraits.hh>
@@ -502,6 +503,13 @@ namespace Fem
///////////////////////////////////////////////////////////////////////////
// DiscreteFunctionSpaceSelector
///////////////////////////////////////////////////////////////////////////
+#ifdef USE_BASEFUNCTIONSET_CODEGEN
+ template <class T>
+ using Storage = Dune::Fem::CodegenStorage<T>;
+#else
+ template <class T>
+ using Storage = Dune::Fem::CachingStorage<T>;
+#endif
template< class FunctionSpaceImp, class GridPartImp, int polOrder, DiscreteFunctionSpaces::Enum dfType, Galerkin::Enum opType >
struct DiscreteFunctionSpaceSelector;
@@ -509,34 +517,48 @@ namespace Fem
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::lagrange, Galerkin::Enum::cg >
{
- typedef LagrangeDiscreteFunctionSpace< FunctionSpaceImp, GridPartImp, polOrder, CachingStorage > type;
+ typedef LagrangeDiscreteFunctionSpace< FunctionSpaceImp, GridPartImp, polOrder, Storage > type;
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::legendre, Galerkin::Enum::dg >
{
- typedef LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, CachingStorage > type;
+ typedef LegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, Storage > type;
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::orthonormal, Galerkin::Enum::dg >
{
- typedef DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, CachingStorage > type;
+ typedef DiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, Storage > type;
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::lagrange, Galerkin::Enum::dg >
{
- typedef LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, CachingStorage > type;
+ typedef LagrangeDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, Storage > type;
};
template< class FunctionSpaceImp, class GridPartImp, int polOrder >
struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::hierarchic_legendre, Galerkin::Enum::dg >
{
- typedef HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, CachingStorage > type;
+ typedef HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpaceImp, GridPartImp, polOrder, Storage > type;
+ };
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+ template< class FunctionSpaceImp, class GridPartImp, int polOrder >
+ struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::gausslobatto, Galerkin::Enum::dg >
+ {
+ typedef FixedOrderDGLagrangeSpace< FunctionSpaceImp, GridPartImp, polOrder, Dune::GaussLobattoPointSet, Storage > type;
};
+ template< class FunctionSpaceImp, class GridPartImp, int polOrder >
+ struct DiscreteFunctionSpaceSelector< FunctionSpaceImp, GridPartImp, polOrder, DiscreteFunctionSpaces::Enum::gausslegendre, Galerkin::Enum::dg >
+ {
+ typedef FixedOrderDGLagrangeSpace< FunctionSpaceImp, GridPartImp, polOrder, Dune::GaussLegendrePointSet, Storage > type;
+ };
+#endif
+
template< class ModelImp, class DiscreteFunctionSpaceImp,
DiffusionFlux::Enum diffFluxId, Formulation::Enum form >
struct DiffusionFluxSelector;
diff --git a/dune/fem-dg/misc/configurator.hh b/dune/fem-dg/misc/configurator.hh
index f8bdbc553144c7786f3e2a4af9966c38345a2099..71bdf13a7ee58f926b42eb56cac99fdc7a68603a 100644
--- a/dune/fem-dg/misc/configurator.hh
+++ b/dune/fem-dg/misc/configurator.hh
@@ -75,6 +75,10 @@ namespace Fem
DiscreteFunctionSpaces::Enum::lagrange;
#elif DISCRETEFUNCTIONSPACESENUM == 5
DiscreteFunctionSpaces::Enum::padaptive;
+#elif DISCRETEFUNCTIONSPACESENUM == 6
+ DiscreteFunctionSpaces::Enum::gausslobatto;
+#elif DISCRETEFUNCTIONSPACESENUM == 7
+ DiscreteFunctionSpaces::Enum::gausslegendre;
#else // DISCRETEFUNCTIONSPACESENUM == 1
DiscreteFunctionSpaces::Enum::orthonormal;
#endif
diff --git a/dune/fem-dg/misc/simulator.hh b/dune/fem-dg/misc/simulator.hh
index 0ef5faa1ea7d273bfaa1acc00a28670e56f9d9bf..d70fff96aba6c6fd8a992095420893663fee70fe 100644
--- a/dune/fem-dg/misc/simulator.hh
+++ b/dune/fem-dg/misc/simulator.hh
@@ -49,12 +49,12 @@
#include <dune/fem-dg/misc/streams.hh>
#if defined USE_BASEFUNCTIONSET_CODEGEN || defined BASEFUNCTIONSET_CODEGEN_GENERATE
-#include <dune/fem/space/basisfunctionset/default_codegen.hh>
#include <dune/fem/space/basisfunctionset/codegen.hh>
#endif
#ifdef USE_BASEFUNCTIONSET_CODEGEN
#warning "Using autogenerated code"
+#include <autogeneratedcode.hh>
#endif
#include <dune/fem-dg/pass/threadhandle.hh>
@@ -269,7 +269,7 @@ namespace Fem
static void apply( const AlgorithmCreator& algCreator, const int polynomialOrder, const bool computeAnyway )
{
#ifdef BASEFUNCTIONSET_CODEGEN_GENERATE
- Dune::Fem::CodegenInfo :: instance().setFileName( autoFilename( CODEDIM, polOrd ) );
+ Dune::Fem::Codegen::CodegenInfo :: instance().setFileName( autoFilename( CODEDIM, polOrd ) );
#endif
if( computeAnyway || polOrd == polynomialOrder )
{
@@ -290,8 +290,9 @@ namespace Fem
static void run( const AlgorithmCreator& algCreator )
{
#ifdef BASEFUNCTIONSET_CODEGEN_GENERATE
+ using namespace Dune::Fem::Codegen;
#if COUNT_FLOPS
- Dune::Fem::CodegenInfo :: instance().setPath(
+ CodegenInfo :: instance().setPath(
Dune::Fem::Parameter::commonInputPath() + "/autogeneratedcodeflops");
#endif
std::cout << "Generating Code \n";
@@ -306,11 +307,11 @@ namespace Fem
Dune::Fem::ForLoop< SimulatePolOrd, MIN_POLORD, MAX_POLORD > :: apply( algCreator, polOrder, bool(MIN_POLORD == MAX_POLORD) );
#ifdef BASEFUNCTIONSET_CODEGEN_GENERATE
}
- catch (const Dune::Fem::CodegenInfoFinished& ) {}
+ catch (const CodegenInfoFinished& ) {}
std::cerr << "Code for k="<< MAX_POLORD << " generated!! " << std::endl;
- Dune::Fem::CodegenInfo :: instance().dumpInfo();
- Dune::Fem::CodegenInfo :: instance().clear();
+ CodegenInfo :: instance().dumpInfo();
+ CodegenInfo :: instance().clear();
std::remove( autoFilename( CODEDIM, MAX_POLORD ).c_str() );
finalizeCodegen();
#endif
diff --git a/dune/fem-dg/models/modelwrapper.hh b/dune/fem-dg/models/modelwrapper.hh
index aed6d44469d37279f01d2f24fa9cc20bebecf5ef..2cabcf41971ec88485b869ea31d63316a2921c06 100644
--- a/dune/fem-dg/models/modelwrapper.hh
+++ b/dune/fem-dg/models/modelwrapper.hh
@@ -198,7 +198,7 @@ namespace Fem
limitedRange_[ i ] = i;
// if method has been filled then modified will be set differently
- advection_.limitedRange( limitedRange_ );
+ advection().limitedRange( limitedRange_ );
}
#ifdef EULER_WRAPPER_TEST
@@ -210,17 +210,15 @@ namespace Fem
void setTime (const double t)
{
BaseType::setTime(t);
-#if HAVE_DUNE_FEMPY
// update model times (only if time method is available on these models)
//! TODO problem without virtualization advection_.setTime(t);
- //! TODO problem without virtualization diffusion_.setTime(t);
+ //! TODO problem without virtualization diffusion().setTime(t);
::detail::CallSetTime< AdvectionModelType,
::detail::CheckTimeMethod< AdvectionModelType >::value >
- ::setTime( const_cast< AdvectionModelType& > (advection_), t );
+ ::setTime( const_cast< AdvectionModelType& > (advection()), t );
::detail::CallSetTime< DiffusionModelType,
::detail::CheckTimeMethod< DiffusionModelType >::value >
- ::setTime( const_cast< DiffusionModelType& > (diffusion_), t );
-#endif
+ ::setTime( const_cast< DiffusionModelType& > (diffusion()), t );
}
double gamma () const { return problem_.gamma(); }
@@ -229,9 +227,9 @@ namespace Fem
void setEntity( const Entity& entity ) const
{
if( hasAdvection )
- advection_.init( entity );
+ advection().init( entity );
if( hasDiffusion )
- diffusion_.init( entity );
+ diffusion().init( entity );
}
inline bool hasStiffSource() const { return AdditionalType::hasStiffSource; }
@@ -242,8 +240,8 @@ namespace Fem
void obtainBounds( RangeType& globalMin, RangeType& globalMax) const
{
- globalMin = 0;
- globalMax = std::numeric_limits<double>::max();
+ assert( hasAdvection );
+ advection().obtainBounds(globalMin, globalMax);
}
bool isConstant( const RangeType& min, const RangeType& max ) const
@@ -259,7 +257,7 @@ namespace Fem
RangeType & s) const
{
assert( hasDiffusion );
- diffusion_.source( local.quadraturePoint(), u, du, s );
+ diffusion().source( local.quadraturePoint(), u, du, s );
return 0;
}
@@ -271,7 +269,7 @@ namespace Fem
RangeType& s) const
{
assert( hasAdvection );
- advection_.source( local.quadraturePoint(), u, du, s );
+ advection().source( local.quadraturePoint(), u, du, s );
return 0;
}
@@ -282,7 +280,7 @@ namespace Fem
JacobianRangeType& result ) const
{
assert( hasAdvection );
- advection_.diffusiveFlux( local.quadraturePoint(), u, du, result );
+ advection().diffusiveFlux( local.quadraturePoint(), u, du, result );
}
template <class LocalEvaluation>
@@ -292,7 +290,7 @@ namespace Fem
{
if( hasDiffusion )
{
- maxValue = diffusion_.diffusionTimeStep( local.entity(), local.quadraturePoint(), 0.0, u );
+ maxValue = diffusion().diffusionTimeStep( local.entity(), local.quadraturePoint(), 0.0, u );
}
}
@@ -301,7 +299,7 @@ namespace Fem
const T& circumEstimate,
const RangeType& u ) const
{
- return diffusion_.diffusionTimeStep( local.entity(), local.quadraturePoint(), circumEstimate, u );
+ return diffusion().diffusionTimeStep( local.entity(), local.quadraturePoint(), circumEstimate, u );
}
// is not used
@@ -314,7 +312,7 @@ namespace Fem
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);
+ advection().linDiffusiveFlux( u, du, local.quadraturePoint(), u, du, A);
}
template <class LocalEvaluation>
@@ -329,7 +327,7 @@ namespace Fem
RangeType u; // fake return variable
const int id = getBoundaryId( local );
// the following fails since is is called with
- return advection_.hasBoundaryValue(id, time(), local.entity(), local.position(), u, u);
+ return advection().hasBoundaryValue(id, time(), local.entity(), local.position(), u, u);
}
// return uRight for insertion into the numerical flux
@@ -342,7 +340,7 @@ namespace Fem
#ifndef NDEBUG
const bool isDirichlet =
#endif
- advection_.boundaryValue(id, time(), local.entity(), local.quadraturePoint(),
+ advection().boundaryValue(id, time(), local.entity(), local.quadraturePoint(),
local.intersection().unitOuterNormal( local.localPosition() ),
uLeft, uRight);
assert( isDirichlet );
@@ -363,7 +361,7 @@ namespace Fem
#ifndef NDEBUG
const bool isFluxBnd =
#endif
- advection_.boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
+ advection().boundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, gLeft);
gLeft *= len;
assert( isFluxBnd );
return 0; // TODO: do something better here
@@ -376,7 +374,7 @@ namespace Fem
JacobianRangeType& diff ) const
{
assert( hasDiffusion );
- diffusion_.diffusiveFlux( local.quadraturePoint(), u, du, diff);
+ diffusion().diffusiveFlux( local.quadraturePoint(), u, du, diff);
}
@@ -395,7 +393,7 @@ namespace Fem
#ifndef NDEBUG
const bool isFluxBnd =
#endif
- diffusion_.diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
+ diffusion().diffusionBoundaryFlux(id, time(), local.entity(), local.quadraturePoint(), normal, uLeft, jacLeft, gLeft);
gLeft *= len;
assert( isFluxBnd );
return 0; // QUESTION: do something better here? Yes, return time step restriction if possible
@@ -411,7 +409,7 @@ namespace Fem
// TODO: add a max speed for the diffusion time step control
// this needs to be added in diffusionTimeStep
assert( hasAdvection );
- advspeed = advection_.maxSpeed( time(), local.entity(), local.quadraturePoint(), unitNormal, u );
+ advspeed = advection().maxSpeed( time(), local.entity(), local.quadraturePoint(), unitNormal, u );
totalspeed = advspeed;
}
@@ -425,7 +423,7 @@ namespace Fem
inline DomainType velocity (const LocalEvaluation& local,
RangeType& u) const
{
- return advection_.velocity( time(), local.entity(), local.quadraturePoint(), u );
+ return advection().velocity( time(), local.entity(), local.quadraturePoint(), u );
}
/////////////////////////////////////////////////////////////////
@@ -437,7 +435,7 @@ namespace Fem
const RangeType& u,
DomainType& velocity) const
{
- velocity = advection_.velocity( time(), en, x, u );
+ velocity = advection().velocity( time(), en, x, u );
}
// we have physical check for this model
@@ -452,7 +450,7 @@ namespace Fem
const DomainType& x,
const RangeType& u) const
{
- return advection_.physical( entity, x, u ) > 0;
+ return advection().physical( entity, x, u ) > 0;
}
// adjust average value if necessary
@@ -463,7 +461,7 @@ namespace Fem
RangeType& u ) const
{
// nothing to be done here for this test case
- advection_.adjustAverageValue( entity, xLocal, u );
+ advection().adjustAverageValue( entity, xLocal, u );
}
// calculate jump between left and right value
@@ -474,7 +472,7 @@ namespace Fem
const RangeType& uRight,
RangeType& jump) const
{
- jump = advection_.jump( it, x, uLeft, uRight );
+ jump = advection().jump( it, x, uLeft, uRight );
}
// calculate jump between left and right value
@@ -485,7 +483,7 @@ namespace Fem
const RangeType& uRight,
RangeType& indicator) const
{
- indicator = advection_.jump( it, x, uLeft, uRight );
+ indicator = advection().jump( it, x, uLeft, uRight );
}
// misc for eoc calculation
@@ -495,8 +493,13 @@ namespace Fem
}
protected:
- const AdvectionModelType& advection_;
- const DiffusionModelType& diffusion_;
+ const AdvectionModelType& advection () const { return advection_; }
+ const DiffusionModelType& diffusion () const { return diffusion_; }
+
+ // store a copy of the models here for thread parallel runs
+ // where we need class variables to be thread safe
+ AdvectionModelType advection_;
+ DiffusionModelType diffusion_;
const ProblemType& problem_;
LimitedRangeType limitedRange_;
diff --git a/dune/fem-dg/operator/dg/CMakeLists.txt b/dune/fem-dg/operator/dg/CMakeLists.txt
index 32a1bc6768f3be52fa16629000f1e3e381d20171..504bdc9b28099bdd78985769e4b269981de0982d 100644
--- a/dune/fem-dg/operator/dg/CMakeLists.txt
+++ b/dune/fem-dg/operator/dg/CMakeLists.txt
@@ -1 +1 @@
-dune_install( discretemodelcommon.hh fluxdiscretemodel.hh fluxoperator.hh operatorbase.hh operatortraits.hh primaldiscretemodel.hh primaloperator.hh)
+dune_install( defaultquadrature.hh discretemodelcommon.hh fluxdiscretemodel.hh fluxoperator.hh operatorbase.hh operatortraits.hh primaldiscretemodel.hh primaloperator.hh)
diff --git a/dune/fem-dg/operator/dg/defaultquadrature.hh b/dune/fem-dg/operator/dg/defaultquadrature.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dd73d6af3ad848632dcdaaec05c69ef22b6a9b63
--- /dev/null
+++ b/dune/fem-dg/operator/dg/defaultquadrature.hh
@@ -0,0 +1,83 @@
+#ifndef DUNE_FEMDG_DEFAULTQUADRATURE_HH
+#define DUNE_FEMDG_DEFAULTQUADRATURE_HH
+
+#include <dune/fem/space/lagrange.hh>
+#include <dune/fem/quadrature/defaultquadratures.hh>
+#include <dune/fem/quadrature/interpolationquadrature.hh>
+
+namespace Dune
+{
+ namespace Fem
+ {
+
+ struct DefaultQuadratureGauss
+ {
+ template <class F, int d>
+ using DefaultQuadratureTraits = Dune::Fem::DefaultQuadratureTraits< F, d >;
+
+ static int volumeOrder( const int polOrder ) { return 2 * polOrder; }
+ static int faceOrder( const int polOrder ) { return 2 * polOrder + 1; }
+ };
+
+ template < class Space >
+ struct DefaultQuadrature : public DefaultQuadratureGauss
+ {
+ };
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+ struct DefaultQuadratureGaussLobatto
+ {
+ template <class F, int d>
+ using DefaultQuadratureTraits = Dune::Fem::GaussLobattoQuadratureTraits< F, d >;
+
+ static int volumeOrder( const int polOrder ) { return (polOrder > 0) ? (2 * polOrder - 1) : 0; }
+ static int faceOrder( const int polOrder ) { return (polOrder > 0) ? (2 * polOrder - 1) : 0; }
+ };
+
+ template < class FunctionSpace, class GridPart, unsigned int order, template< class > class Storage >
+ struct DefaultQuadrature< Dune::Fem::FixedOrderDGLagrangeSpace< FunctionSpace, GridPart, order, Dune::GaussLobattoPointSet, Storage > >
+ : public DefaultQuadratureGaussLobatto
+ {
+ };
+
+ template< class LFEMap >
+ struct DefaultQuadratureSpec : DefaultQuadratureGauss
+ {};
+
+ template< class FunctionSpace, class GridPart, unsigned int order >
+ struct DefaultQuadratureSpec< FixedOrderLagrangeFiniteElementMap< FunctionSpace, GridPart, order, Dune::GaussLobattoPointSet > >
+ : public DefaultQuadratureGaussLobatto
+ {
+ };
+
+ struct DefaultQuadratureGaussLegendre
+ {
+ template <class F, int d>
+ using DefaultQuadratureTraits = Dune::Fem::GaussLegendreQuadratureTraits< F, d >;
+
+ static int volumeOrder( const int polOrder ) { return 2 * polOrder + 1; }
+ static int faceOrder( const int polOrder ) { return 2 * polOrder + 1; }
+ };
+
+ template < class FunctionSpace, class GridPart, unsigned int order, template< class > class Storage >
+ struct DefaultQuadrature< Dune::Fem::FixedOrderDGLagrangeSpace< FunctionSpace, GridPart, order, Dune::GaussLegendrePointSet, Storage > >
+ : public DefaultQuadratureGaussLegendre
+ {};
+
+ template< class FunctionSpace, class GridPart, unsigned int order >
+ struct DefaultQuadratureSpec< FixedOrderLagrangeFiniteElementMap< FunctionSpace, GridPart, order, Dune::GaussLegendrePointSet > >
+ : public DefaultQuadratureGaussLegendre
+ {
+ };
+
+ template< class LFEMap, class FunctionSpace, template< class > class Storage >
+ struct DefaultQuadrature< DiscontinuousLocalFiniteElementSpace< LFEMap, FunctionSpace, Storage > >
+ : public DefaultQuadratureSpec< LFEMap >
+ {
+ };
+
+#endif
+ }
+}
+
+#endif
diff --git a/dune/fem-dg/operator/dg/dgpyoperator.hh b/dune/fem-dg/operator/dg/dgpyoperator.hh
index 661ecd7d045393f3a2ff0b8679bb64346f23acea..7239d322390f4e9cb23946ecab812aaf288c5c30 100644
--- a/dune/fem-dg/operator/dg/dgpyoperator.hh
+++ b/dune/fem-dg/operator/dg/dgpyoperator.hh
@@ -8,6 +8,10 @@
#include <dune/fem/solver/timeprovider.hh>
#include <dune/fem/operator/common/spaceoperatorif.hh>
+#include <dune/fem/space/lagrange.hh>
+#include <dune/fem/quadrature/interpolationquadrature.hh>
+
+
#include <dune/fem-dg/algorithm/evolution.hh>
// dune-fem-dg includes
@@ -18,10 +22,10 @@
#include <dune/fem-dg/models/modelwrapper.hh>
#include <dune/fem-dg/misc/algorithmcreatorselector.hh>
#include <dune/fem-dg/operator/adaptation/estimator.hh>
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
-#if HAVE_DUNE_FEMPY
+// this is part of dune-fem now
#include <dune/fempy/quadrature/fempyquadratures.hh>
-#endif
namespace Dune
{
@@ -79,14 +83,16 @@ namespace Fem
typedef typename DiffusionFluxSelector< ModelType, DiscreteFunctionSpaceType, diffFluxId, formId >::type DiffusionFluxType;
- typedef DefaultOperatorTraits< ModelType, DestinationType, AdvectionFluxType, DiffusionFluxType,
- std::tuple<>, typename DiscreteFunctionSpaceType::FunctionSpaceType,
-#if HAVE_DUNE_FEMPY
- Dune::FemPy::FempyQuadratureTraits, // use quadratures from dune-fempy
-#else
- Dune::Fem::DefaultQuadratureTraits,
-#endif
- threading > OpTraits;
+ typedef typename std::conditional< Additional::defaultQuadrature,
+ DefaultOperatorTraits< ModelType, DestinationType, AdvectionFluxType, DiffusionFluxType,
+ std::tuple<>, typename DiscreteFunctionSpaceType::FunctionSpaceType,
+ DefaultQuadrature< DiscreteFunctionSpaceType >::template DefaultQuadratureTraits,
+ threading>,
+ // fempy quadratures
+ DefaultOperatorTraits< ModelType, DestinationType, AdvectionFluxType, DiffusionFluxType,
+ std::tuple<>, typename DiscreteFunctionSpaceType::FunctionSpaceType,
+ Dune::FemPy::FempyQuadratureTraits,
+ threading> > ::type OpTraits;
typedef AdvectionDiffusionOperatorSelector< OpTraits, formId, limiterId > OperatorSelectorType ;
@@ -98,6 +104,7 @@ namespace Fem
typedef Estimator< DestinationType, typename ModelType::ProblemType > GradientIndicatorType ;
typedef AdaptIndicator< IndicatorType, GradientIndicatorType > AdaptIndicatorType;
+ typedef typename FullOperatorType :: TroubledCellIndicatorType TroubledCellIndicatorType;
// solver selection, available fem, istl, petsc, ...
typedef typename MatrixFreeSolverSelector< solverId, symmetric > :: template LinearInverseOperatorType< DiscreteFunctionSpaceType, DiscreteFunctionSpaceType > LinearSolverType ;
@@ -114,9 +121,9 @@ namespace Fem
model_( advectionModel, diffusionModel, problem_ ),
advFluxPtr_(),
advFlux_( advectionFlux( parameter, std::integral_constant< bool, fluxIsUserDefined >() )),
- fullOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
- explOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
- implOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
+ fullOperator_( space.gridPart(), model_, advFlux_, extra_, "full", parameter ),
+ explOperator_( space.gridPart(), model_, advFlux_, extra_, "expl", parameter ),
+ implOperator_( space.gridPart(), model_, advFlux_, extra_, "impl", parameter ),
adaptIndicator_( space, model_, advFlux_, extra_, name()+"_adaptind", parameter )
{
}
@@ -132,9 +139,9 @@ namespace Fem
model_( advectionModel, diffusionModel, problem_ ),
advFluxPtr_(),
advFlux_( advFlux ),
- fullOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
- explOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
- implOperator_( space.gridPart(), model_, advFlux_, extra_, name(), parameter ),
+ fullOperator_( space.gridPart(), model_, advFlux_, extra_, "full", parameter ),
+ explOperator_( space.gridPart(), model_, advFlux_, extra_, "expl", parameter ),
+ implOperator_( space.gridPart(), model_, advFlux_, extra_, "impl", parameter ),
adaptIndicator_( space, model_, advFlux_, extra_, name()+"_adaptind", parameter )
{
}
@@ -183,6 +190,13 @@ namespace Fem
}
}
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator)
+ {
+ fullOperator_.setTroubledCellIndicator(indicator);
+ explOperator_.setTroubledCellIndicator(indicator);
+ }
+
+
/** \copydoc SpaceOperatorInterface::setTime */
void setTime( const double time )
{
@@ -191,6 +205,23 @@ namespace Fem
double timeStepEstimate() const { return fullOperator_.timeStepEstimate(); }
+ //! return number of interior elements visited by the operator
+ inline size_t gridSizeInterior() const
+ {
+ size_t elem = std::max( fullOperator_.numberOfElements(),
+ explOperator_.numberOfElements() );
+
+ // this can occur if element size if requested
+ // before operator had been applied
+ if( elem == 0 )
+ {
+ const auto end = space_.end();
+ size_t elem = 0;
+ for( auto it = space_.begin(); it != end; ++it, ++elem );
+ }
+ return elem;
+ }
+
//// End Methods from SpaceOperatorInterface /////
protected:
@@ -222,7 +253,6 @@ namespace Fem
mutable double fixedTimeStep_ ;
mutable bool initialized_;
-
};
} // end namespace Fem
diff --git a/dune/fem-dg/operator/dg/discretemodelcommon.hh b/dune/fem-dg/operator/dg/discretemodelcommon.hh
index cdf1686e2fff36189940333fd02e41adcc8ec058..92006147073e8dc6dbc75b94b935accc5366636f 100644
--- a/dune/fem-dg/operator/dg/discretemodelcommon.hh
+++ b/dune/fem-dg/operator/dg/discretemodelcommon.hh
@@ -304,7 +304,11 @@ namespace Fem
return hasBndValue;
}
- const ModelType& model_;
+ // store an instance here so that for thread parallel
+ // runs class variables are thread private since discrete models
+ // are thread private
+ ModelType model_;
+
const AdvectionFluxType& numflux_;
mutable RangeType uBnd_;
}; /*@LST0E@*/
diff --git a/dune/fem-dg/operator/dg/operatorbase.hh b/dune/fem-dg/operator/dg/operatorbase.hh
index 1a9521a7dfa5f169ccdc3f0dfd07dcf255204f59..74bef4e86e2ae64098cf3dda5ac5f80653ea9ceb 100644
--- a/dune/fem-dg/operator/dg/operatorbase.hh
+++ b/dune/fem-dg/operator/dg/operatorbase.hh
@@ -92,6 +92,11 @@ namespace Fem
public:
using BaseType::operator () ;
+ // dummy method for a troubled cell indicator to be passed to the
+ // limited advection operator
+ typedef void* TroubledCellIndicatorType;
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator) {}
+
typedef typename Traits::ModelType ModelType;
typedef typename ModelType::ProblemType ProblemType ;
diff --git a/dune/fem-dg/operator/dg/operatortraits.hh b/dune/fem-dg/operator/dg/operatortraits.hh
index 55b21a764661416243bda9653c1f6f7df2dc2c52..23f4ac9a2cfd4f12e4d554ecb4ffe00f0c9b5363 100644
--- a/dune/fem-dg/operator/dg/operatortraits.hh
+++ b/dune/fem-dg/operator/dg/operatortraits.hh
@@ -1,56 +1,29 @@
#ifndef DUNE_FEM_DG_OPERATORTRAITS_HH
#define DUNE_FEM_DG_OPERATORTRAITS_HH
+#include <dune/fem/quadrature/dunequadratures.hh>
#include <dune/fem/quadrature/cachingquadrature.hh>
+#include <dune/fem/quadrature/interpolationquadrature.hh>
#include <dune/fem-dg/operator/fluxes/diffusion/fluxes.hh>
#include <dune/fem-dg/operator/adaptation/adaptation.hh>
#include <dune/fem/space/finitevolume/space.hh>
#include <dune/fem/function/adaptivefunction/adaptivefunction.hh>
+
#include <dune/fem-dg/operator/fluxes/diffusion/parameters.hh>
+
+#if HAVE_DUNE_POLYGONGRID
+#include <dune/polygongrid/declaration.hh>
+#endif
+
+#if HAVE_OPM_GRID
+#include <opm/grid/polyhedralgrid/declaration.hh>
+#endif // #if HAVE_OPM_GRID
+
namespace Dune
{
namespace Fem
{
-
- // traits for the operator passes
- template< class GridPart,
- int polOrd,
- class ModelImp,
- class DiscreteFunctionImp,
- class AdvectionFluxImp,
- class DiffusionFluxImp,
- class LimiterIndicatorFunctionImp,
- class AdaptationHandlerImp,
- class ExtraParameterTupleImp = std::tuple<>,
- template <class F, int d> class QuadratureTraits = Dune::Fem::DefaultQuadratureTraits
- >
- struct OperatorTraits
- {
- typedef GridPart GridPartType;
- typedef typename GridPartType::GridType GridType;
-
- typedef ModelImp ModelType;
- typedef AdvectionFluxImp AdvectionFluxType;
- typedef DiffusionFluxImp DiffusionFluxType;
-
- // polynomial order of ansatz space
- static const int polynomialOrder = polOrd;
-
- typedef DiscreteFunctionImp DestinationType ;
- //static_assert( std::is_same<typename ModelType::RangeType, typename DiscreteFunctionType::RangeType>::value, "range type does not fit.");
- typedef typename DestinationType::DiscreteFunctionSpaceType DiscreteFunctionSpaceType;
-
- typedef Fem::CachingQuadrature< GridPartType, 0, QuadratureTraits > VolumeQuadratureType;
- typedef Fem::CachingQuadrature< GridPartType, 1, QuadratureTraits > FaceQuadratureType;
-
- typedef LimiterIndicatorFunctionImp LimiterIndicatorType;
-
- typedef AdaptationHandlerImp AdaptationHandlerType ;
-
- typedef ExtraParameterTupleImp ExtraParameterTupleType;
- };
-
// traits for the operator passes
template< class ModelImp,
class DiscreteFunctionImp,
@@ -58,7 +31,7 @@ namespace Fem
class DiffusionFluxImp,
class ExtraParameterTupleImp = std::tuple<>,
class AdaptationHandlerFunctionSpaceImp = typename DiscreteFunctionImp::DiscreteFunctionSpaceType::FunctionSpaceType,
- template <class F, int d> class QuadratureTraits = Dune::Fem::DefaultQuadratureTraits,
+ template <class F, int d> class QuadratureTraits = DefaultQuadrature< typename DiscreteFunctionImp::DiscreteFunctionSpaceType >::template DefaultQuadratureTraits,
bool enableThreaded =
// static cmake variables provided by dune-fem
#ifdef USE_SMP_PARALLEL
@@ -89,9 +62,8 @@ namespace Fem
static_assert( std::is_same<typename ModelType::RangeType, typename DestinationType::RangeType>::value, "range type does not fit.");
-
// default quadrature selection should be CachingQuadrature
- template <class DFS>
+ template <class Grid, bool caching>
struct SelectQuadrature
{
typedef Fem::CachingQuadrature< GridPartType, 0, QuadratureTraits > VolumeQuadratureType;
@@ -100,15 +72,46 @@ namespace Fem
// if selected discrete function space has no caching storage,
// then select ElementQuadrature, and not CachingQuadrature
- template <class FS, class GP, int polOrd,
- template <class, class, int, template <class> class> class DFS>
- struct SelectQuadrature< DFS< FS, GP, polOrd, Dune::Fem::SimpleStorage > >
+ // also, if the grid is polygonal or polyhedral
+ // CachingQuadarature cannot be used
+ template <class Grid>
+ struct SelectQuadrature< Grid, false >
{
typedef Fem::ElementQuadrature< GridPartType, 0, QuadratureTraits > VolumeQuadratureType;
typedef Fem::ElementQuadrature< GridPartType, 1, QuadratureTraits > FaceQuadratureType;
};
- typedef SelectQuadrature< DiscreteFunctionSpaceType > SelectQuadratureType;
+ template <class Grid>
+ struct CheckGrid { static const bool value = true; };
+
+ // also, if the grid is polygonal or polyhedral
+ // CachingQuadarature cannot be used
+#if HAVE_DUNE_POLYGONGRID
+ template <class ct>
+ struct CheckGrid< Dune::PolygonGrid< ct > > { static const bool value = false; };
+#endif
+
+#if HAVE_OPM_GRID
+ template <int dim, int dimworld, class ct>
+ struct CheckGrid< Dune::PolyhedralGrid< dim, dimworld, ct > > { static const bool value = false; };
+#endif
+
+ template <class DFS>
+ struct CheckSpace { static const bool value = true; };
+
+ // if selected discrete function space has no caching storage,
+ // then select ElementQuadrature, and not CachingQuadrature
+ template <class FS, class GP, int polOrd,
+ template <class, class, int, template <class> class> class DFS>
+ struct CheckSpace< DFS< FS, GP, polOrd, Dune::Fem::SimpleStorage > >
+ {
+ static const bool value = false;
+ };
+
+ // define quad selector which contains the correct types for quadratures
+ typedef SelectQuadrature< GridType,
+ CheckSpace< DiscreteFunctionSpaceType > :: value &&
+ CheckGrid < GridType > :: value > SelectQuadratureType;
typedef typename SelectQuadratureType::VolumeQuadratureType VolumeQuadratureType;
typedef typename SelectQuadratureType::FaceQuadratureType FaceQuadratureType;
diff --git a/dune/fem-dg/operator/dg/primaloperator.hh b/dune/fem-dg/operator/dg/primaloperator.hh
index 8bd6882b824d5beebe45ac481bd40a1ef5e5470e..8946d663ec9b7344b2aab3627f21d1bfc9c47e79 100644
--- a/dune/fem-dg/operator/dg/primaloperator.hh
+++ b/dune/fem-dg/operator/dg/primaloperator.hh
@@ -320,6 +320,7 @@ namespace Fem
// Destination of Limiter can differ from the operators destination type depending on FV or DG
// standard limiter and scaling limiter
typedef Limiter< DestinationType, LimiterDestinationType, LimiterDiscreteModelType, threading, ModelType::scalingLimiter > LimiterType;
+ typedef typename LimiterType :: TroubledCellIndicatorType TroubledCellIndicatorType;
// select non-blocking communication handle
typedef typename
@@ -382,17 +383,20 @@ namespace Fem
, indicator_()
, diffFlux_( gridPart_, model_, DGPrimalDiffusionFluxParameters( ParameterKey::generate( name, "dgdiffusionflux." ), parameter ) )
, discreteModel1_( model_, advFlux_, diffFlux_ )
- , limiter_( space_, limiterSpace_, model_ )
+ , limiter_( space_, limiterSpace_, model_, parameter )
, pass0_()
, pass2_( discreteModel1_, pass0_, space_ )
, counter_(0)
, limitTime_( 0 )
, computeTime_( 0 )
+ , name_( name )
+ , verbose_( Dune::Fem::Parameter::verbose() ) // ||parameter.verbose() )
{
}
virtual ~DGLimitedAdvectionOperator() {
- std::cout << "~DGLimitedAdvectionOperator: op calls = " << counter_ << " T_l = " << limitTime_ << " T_op = " << computeTime_ << std::endl;
+ if( verbose_ )
+ std::cout << "~DGLimitedAdvectionOperator("<<name_<<"): op calls = " << counter_ << " T_l = " << limitTime_ << " T_op = " << computeTime_ << std::endl;
}
void activateLinear() const {
@@ -460,7 +464,7 @@ namespace Fem
inline size_t numberOfElements () const
{
- return pass2_.numberOfElements();
+ return std::max( pass2_.numberOfElements(), limiter_.numberOfElements() );
}
/*
@@ -481,6 +485,11 @@ namespace Fem
LimiterCall< LimiterType, polOrd >::limit( limiter_, arg, U );
}
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator)
+ {
+ limiter_.setTroubledCellIndicator(indicator);
+ }
+
void printmyInfo(std::string filename) const
{
std::ostringstream filestream;
@@ -513,11 +522,11 @@ namespace Fem
}
protected:
- GridPartType& gridPart_;
- const ModelType& model_;
- const AdvectionFluxType& advFlux_;
+ GridPartType& gridPart_;
+ const ModelType& model_;
+ const AdvectionFluxType& advFlux_;
- SpaceType space_;
+ SpaceType space_;
LimiterSpaceType limiterSpace_;
mutable LimiterDestinationType limitedU_;
@@ -539,6 +548,9 @@ namespace Fem
mutable double limitTime_;
mutable double computeTime_;
+
+ std::string name_;
+ const bool verbose_;
};
diff --git a/dune/fem-dg/operator/fluxes/advection/llfadvflux.hh b/dune/fem-dg/operator/fluxes/advection/llfadvflux.hh
index 593ac493f3d775a14eb9a75ae36de79a8a4e1e0e..f7326d377c5c1db5070a92035901b840eacf1391 100644
--- a/dune/fem-dg/operator/fluxes/advection/llfadvflux.hh
+++ b/dune/fem-dg/operator/fluxes/advection/llfadvflux.hh
@@ -70,20 +70,26 @@ namespace Fem
const auto faceArea = normal.two_norm();
normal *= 1./faceArea;
+ double maxspeedl, maxspeedr;
+ double viscparal, viscparar;
+
FluxRangeType anaflux;
+ // assume that model is set to left entity
+ model_.setEntity( left.entity() );
+
model_.advection( left, uLeft, jacLeft, anaflux );
// set gLeft
anaflux.mv( normal, gLeft );
+ model_.maxSpeed( left, normal, uLeft, viscparal, maxspeedl );
+
+ model_.setEntity( right.entity() );
+
model_.advection( right, uRight, jacRight, anaflux );
// add to gLeft
anaflux.umv( normal, gLeft );
- double maxspeedl, maxspeedr;
- double viscparal, viscparar;
-
- model_.maxSpeed( left, normal, uLeft, viscparal, maxspeedl );
model_.maxSpeed( right, normal, uRight, viscparar, maxspeedr );
const double maxspeed = std::max( maxspeedl, maxspeedr);
@@ -99,6 +105,8 @@ namespace Fem
// conservation property
gRight = gLeft;
+ model_.setEntity( left.entity() );
+
return maxspeed * faceArea;
}
};
diff --git a/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh b/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
index 764fd1bc05798bc7028a9acdc92f6ba779293b54..f9c7e4be3c626a4c4c0011889c1f6f6c3e04c154 100644
--- a/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
+++ b/dune/fem-dg/operator/fluxes/diffusion/dgprimalfluxes.hh
@@ -11,6 +11,8 @@
#include <dune/fem/storage/dynamicarray.hh>
#include <dune/fem-dg/pass/dgmasspass.hh>
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
+
#include "fluxbase.hh"
namespace Dune
@@ -62,7 +64,8 @@ namespace Fem
typedef typename DiscreteGradientSpaceType::RangeType GradientType;
typedef Fem::TemporaryLocalFunction< DiscreteGradientSpaceType > LiftingFunctionType;
- typedef Fem::CachingQuadrature< GridPartType, 0> VolumeQuadratureType ;
+ typedef Fem::CachingQuadrature< GridPartType, 0,
+ DefaultQuadrature<DiscreteGradientSpaceType>::template DefaultQuadratureTraits> VolumeQuadratureType ;
typedef Fem::LocalMassMatrix
< DiscreteGradientSpaceType, VolumeQuadratureType > LocalMassMatrixType;
@@ -96,7 +99,7 @@ namespace Fem
#ifdef USE_CACHED_INVERSE_MASSMATRIX
, localMassMatrix_( InverseMassProviderType :: getObject( KeyType( gradSpc_.gridPart() ) ) )
#else
- , localMassMatrix_( gradSpc_, 2*gradSpc_.order() )
+ , localMassMatrix_( gradSpc_, [this](const int order) { return DefaultQuadrature<DiscreteGradientSpaceType>::volumeOrder(order); } )
#endif
, isInitialized_( 0 )
{
diff --git a/dune/fem-dg/operator/limiter/indicatorbase.hh b/dune/fem-dg/operator/limiter/indicatorbase.hh
new file mode 100644
index 0000000000000000000000000000000000000000..10edaf26ce387cf429454cb1a5ea2dbffdc4ae0a
--- /dev/null
+++ b/dune/fem-dg/operator/limiter/indicatorbase.hh
@@ -0,0 +1,13 @@
+#pragma once
+namespace Dune
+{
+ namespace Fem
+ {
+ template< class DiscreteFunction >
+ struct TroubledCellIndicatorBase
+ {
+ typedef typename DiscreteFunction::LocalFunctionType LocalFunctionType;
+ virtual double operator()( const DiscreteFunction& U, const LocalFunctionType& uEn ) const = 0;
+ };
+ }
+}
diff --git a/dune/fem-dg/operator/limiter/limiter.hh b/dune/fem-dg/operator/limiter/limiter.hh
index 5dafcdb848d263c6282abcc0e8f4bd4a663081ab..804e089f180a37385a0449b78f179f6fae7a2e92 100644
--- a/dune/fem-dg/operator/limiter/limiter.hh
+++ b/dune/fem-dg/operator/limiter/limiter.hh
@@ -115,9 +115,12 @@ namespace Fem
ThreadPass < InnerPassType, ThreadIteratorType, true>,
InnerPassType > :: type LimitPassType;
- typedef typename LimiterDiscreteModelType::IndicatorType IndicatorType;
+ typedef typename LimiterDiscreteModelType :: IndicatorType IndicatorType;
typedef typename IndicatorType :: DiscreteFunctionSpaceType IndicatorSpaceType;
+ // InnerPassType == LimitPass in case of threading
+ typedef typename InnerPassType :: TroubledCellIndicatorType TroubledCellIndicatorType;
+
public:
Limiter( const DomainSpaceType& domainSpace,
const double lowerBound, const double upperBound )
@@ -133,15 +136,16 @@ namespace Fem
Limiter( const DomainSpaceType& domainSpace,
const RangeSpaceType& rangeSpace,
- const ModelType& model )
+ const ModelType& model,
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container() )
: domainSpace_( domainSpace )
, rangeSpace_( rangeSpace )
, indicatorSpace_( domainSpace_.gridPart() )
, indicator_("indicator", indicatorSpace_ )
, model_( model )
- , discreteModel_( model_, domainSpace_.order() )
+ , discreteModel_( model_, domainSpace_.order(), parameter )
, startPass_()
- , limitPass_( discreteModel_ , startPass_, rangeSpace_ )
+ , limitPass_( discreteModel_ , startPass_, rangeSpace_, parameter )
{
discreteModel_.setIndicator( &indicator_ );
@@ -159,6 +163,11 @@ namespace Fem
limitPass_( arg, dest );
}
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator)
+ {
+ limitPass_.setTroubledCellIndicator(indicator);
+ }
+
bool activated( const EntityType& entity ) const
{
return std::abs( indicator_.localFunction( entity )[ 0 ] ) > 0.0;
@@ -173,6 +182,11 @@ namespace Fem
return limitPass_.computeTimeSteps();
}
+ inline size_t numberOfElements () const
+ {
+ return limitPass_.numberOfElements();
+ }
+
private:
Limiter( const Limiter& );
diff --git a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
index 2601492793459f8a753939ea42b9fb1c498f9be7..a8796ff78f8e03511340243af2bfcec504538e8e 100644
--- a/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
+++ b/dune/fem-dg/operator/limiter/limiterdiscretemodel.hh
@@ -439,12 +439,13 @@ namespace Fem
//! returns true, if we have an inflow boundary
template <class ArgumentTuple>
bool checkDirection(const IntersectionType& it,
- const double time, const FaceLocalDomainType& x,
+ const FaceLocalDomainType& xFace,
+ const LocalDomainType& xLocal,
const ArgumentTuple& uLeft) const
{
// evaluate velocity
- model_.velocity(this->inside(), this->inside().geometry().local( it.geometry().global(x) ), uLeft[ uVar ], velocity_);
- return checkDirection(it, x, velocity_);
+ model_.velocity(this->inside(), xLocal, uLeft[ uVar ], velocity_);
+ return checkDirection(it, xFace, velocity_);
}
protected:
diff --git a/dune/fem-dg/operator/limiter/limiterutility.hh b/dune/fem-dg/operator/limiter/limiterutility.hh
index 899c33368a9c6e70b76c1fb3e0b805b99bc31a75..1043f1c155c9d6d2c4def8ccd182dce7e17828f3 100644
--- a/dune/fem-dg/operator/limiter/limiterutility.hh
+++ b/dune/fem-dg/operator/limiter/limiterutility.hh
@@ -287,13 +287,14 @@ namespace Fem
}
- template <class GridPart, class Entity, class EvalAverage>
+ template <class GridPart, class Entity, class EvalAverage, class LFEn>
static void
setupNeighborValues( const GridPart& gridPart,
const Entity& entity,
const EvalAverage& average,
const DomainType& entityCenter,
const RangeType& entityValue,
+ const LFEn& lfEn,
const bool StructuredGrid,
Flags& flags,
std::vector< DomainType >& baryCenters,
@@ -320,7 +321,7 @@ namespace Fem
const auto endit = gridPart.iend( entity );
- const bool cartesianGrid = flags.cartesian;
+ //const bool cartesianGrid = flags.cartesian;
// loop over all neighbors
for (auto it = gridPart.ibegin( entity ); it != endit; ++it )
@@ -394,7 +395,12 @@ namespace Fem
point -= entityCenter;
const double length = (point * lambda);
- const double factorLength = (cartesianGrid) ? 2.0 * length : length ;
+ // this version mirrors the point outside of the domain but only
+ // for Cartesian grids - I think having different behaviour
+ // depending on the grid not so satisfactory
+ // const double factorLength = (cartesianGrid) ? 2.0 * length : length ;
+ // instead put the point onto the boundary
+ const double factorLength = length ;
lambda *= factorLength;
assert( lambda.two_norm () > 0 );
@@ -409,7 +415,13 @@ namespace Fem
const DomainType pointOnBoundary = lambda + entityCenter;
// evaluate data on boundary
- if( average.boundaryValue( entity, intersection, interGeo, pointOnBoundary, entityValue, jumpNeighborEntity ) )
+ // This version uses the bary center value to construct the
+ // boundary value
+ // if( average.boundaryValue( entity, intersection, interGeo, pointOnBoundary, entityValue, jumpNeighborEntity ) )
+ // instead evaluate the discrete function directly on the boundary
+ RangeType enBnd(entityValue);
+ lfEn.evaluate(entity.geometry().local(pointOnBoundary), enBnd);
+ if( average.boundaryValue( entity, intersection, interGeo, pointOnBoundary, enBnd, jumpNeighborEntity ) )
{
jumpNeighborEntity -= entityValue;
}
@@ -924,9 +936,9 @@ namespace Fem
protected:
void printInfo(const std::string& name ) const
{
- //if( Parameter::verbose() )
+ if( Parameter::verbose() )
{
- std::cout << "LimiterFunction: " << name << " with epsilon = " << epsilon_ << std::endl;
+ std::cout << "LimiterFunction (" << name << ") with femdg.limiter.epsilon: " << epsilon_ << std::endl;
}
}
};
diff --git a/dune/fem-dg/operator/limiter/limitpass.hh b/dune/fem-dg/operator/limiter/limitpass.hh
index cf1ee086d2a557348177c20a22f4c59d595ec0c4..5605e0e8f3b81c43946a7151fd6adfacdd7e8c2e 100644
--- a/dune/fem-dg/operator/limiter/limitpass.hh
+++ b/dune/fem-dg/operator/limiter/limitpass.hh
@@ -13,8 +13,6 @@
#include <dune/fem/quadrature/cornerpointset.hh>
#include <dune/fem/io/parameter.hh>
-#include <dune/fem/operator/1order/localmassmatrix.hh>
-
#include <dune/fem/space/common/adaptationmanager.hh>
#include <dune/fem/space/common/basesetlocalkeystorage.hh>
#include <dune/fem/space/common/capabilities.hh>
@@ -23,6 +21,7 @@
#include <dune/fem/space/finitevolume.hh>
#include <dune/fem/function/adaptivefunction.hh>
+#include <dune/fem/function/localfunction/bindable.hh>
#include <dune/fem/misc/compatibility.hh>
#include <dune/fem/misc/threads/threadmanager.hh>
@@ -36,10 +35,10 @@
#include <dune/fem-dg/operator/limiter/limiterdiscretemodel.hh>
#include <dune/fem-dg/operator/limiter/lpreconstruction.hh>
-#if HAVE_DUNE_OPTIM
-#define WANT_DUNE_OPTIM 1
-//#define WANT_DUNE_OPTIM 0
-#endif
+#include <dune/fem-dg/operator/limiter/smoothness.hh>
+#include <dune/fem-dg/operator/limiter/indicatorbase.hh>
+
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
//*************************************************************
@@ -88,7 +87,8 @@ namespace Fem
// call problem checkDirection
const typename BaseType::EntityType &entity = intersection.inside();
- return discreteModel().checkPhysical( entity, entity.geometry().local( intersection.geometry().global( quadrature.localPoint( qp ) ) ), ranges_ );
+ //return discreteModel().checkPhysical( entity, entity.geometry().local( intersection.geometry().global( quadrature.localPoint( qp ) ) ), ranges_ );
+ return discreteModel().checkPhysical( entity, quadrature.point( qp ), ranges_ );
}
// check whether we have inflow or outflow direction
@@ -102,7 +102,7 @@ namespace Fem
assert( quadPoint_ == qp );
// call checkDirection() on discrete model
- return discreteModel().checkDirection( intersection, time(), quadrature.localPoint( qp ), ranges_ );
+ return discreteModel().checkDirection( intersection, quadrature.localPoint( qp ), quadrature.point( qp ), ranges_ );
}
protected:
using BaseType::discreteModel;
@@ -253,7 +253,8 @@ namespace Fem
typedef AdaptationMethod<GridType> AdaptationMethodType;
//! id for choosing admissible linear functions
- enum AdmissibleFunctions { DGFunctions = 0, ReconstructedFunctions = 1 , BothFunctions = 2 };
+ //! _default refers to muscl(1) on simplices and LP limiter (3) on cubes and polygons
+ enum AdmissibleFunctions { dgonly = 0, muscl = 1 , muscldg = 2, lp = 3, _default = 4 };
//! returns true if pass is currently active in the pass tree
using BaseType :: active ;
@@ -261,8 +262,10 @@ namespace Fem
//! type of Cartesian grid checker
typedef CheckCartesian< GridPartType > CheckCartesianType;
+ //! function describing an external troubled cell indicator
+ typedef std::shared_ptr< TroubledCellIndicatorBase<ArgumentFunctionType> > TroubledCellIndicatorType;
+
protected:
-#if WANT_DUNE_OPTIM
typedef typename GridPartType :: GridViewType GridViewType ;
struct BoundaryValue
@@ -280,7 +283,78 @@ namespace Fem
};
typedef Dune::FV::LPReconstruction< GridViewType, RangeType, BoundaryValue > LinearProgramming;
-#endif
+
+ struct ConstantFunction :
+ public Dune::Fem::BindableGridFunction< GridPartType, Dune::Dim<dimRange> >
+ {
+ typedef Dune::Fem::BindableGridFunction<GridPartType, Dune::Dim<dimRange> > Base;
+ typedef typename Base::RangeType RangeType;
+ using Base::Base;
+
+ RangeType value_;
+
+ template <class Quadrature, class RangeArray>
+ void evaluateQuadrature(const Quadrature& quadrature, RangeArray &values) const
+ {
+ const int nop = quadrature.nop();
+ values.resize( nop );
+ std::fill( values.begin(), values.end(), value_ );
+ }
+
+ template <class Point>
+ void evaluate(const Point &x, typename Base::RangeType &ret) const
+ {
+ ret = value_;
+ }
+
+ bool isConstant() const { return true; }
+
+ unsigned int order() const { return 0; }
+ std::string name() const { return "LimmitPass::ConstantFunction"; } // needed for output
+ };
+
+ struct LinearFunction : public ConstantFunction
+ {
+ DomainType center_;
+ GradientType gradient_;
+
+ typedef ConstantFunction Base;
+ typedef typename Base::RangeType RangeType;
+ using Base::Base;
+ using Base::value_;
+
+ template <class Quadrature, class RangeArray>
+ void evaluateQuadrature(const Quadrature& quadrature, RangeArray &values) const
+ {
+ const int nop = quadrature.nop();
+ values.resize( nop );
+ for( int qp = 0; qp<nop; ++qp )
+ {
+ evaluate( quadrature[ qp ], values[ qp ] );
+ }
+ }
+
+ template <class Point>
+ void evaluate(const Point &x, RangeType &ret) const
+ {
+ Base::evaluate( x, ret );
+
+ // get global coordinate
+ auto point = Base::global( x );
+ point -= center_;
+
+ // evaluate linear function
+ for( int r = 0; r < dimRange; ++r )
+ {
+ ret[ r ] += (gradient_[ r ] * point);
+ }
+ }
+
+ bool isConstant() const { return gradient_.infinity_norm() < 1e-10; }
+
+ unsigned int order() const { return 1; }
+ std::string name() const { return "LimmitPass::LinearFunction"; } // needed for output
+ };
public:
@@ -297,20 +371,39 @@ namespace Fem
PreviousPassType& pass,
const DiscreteFunctionSpaceType& spc,
const int vQ = -1,
- const int fQ = -1 ) :
+ const int fQ = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose() )
+ : LimitDGPass(problem, pass, spc, Dune::Fem::Parameter::container(), vQ, fQ, verbose )
+ {}
+
+ //- Public methods
+ /** \brief constructor
+ *
+ * \param problem Actual problem definition (see problem.hh)
+ * \param pass Previous pass
+ * \param spc Space belonging to the discrete function local to this pass
+ * \param parameter Parameter reader for dynamic parameters
+ * \param vQ order of volume quadrature
+ * \param fQ order of face quadrature
+ */
+
+ LimitDGPass(DiscreteModelType& problem,
+ PreviousPassType& pass,
+ const DiscreteFunctionSpaceType& spc,
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container(),
+ const int vQ = -1,
+ const int fQ = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose() ) :
BaseType(pass, spc),
- caller_( 0 ),
+ caller_(),
discreteModel_(problem),
currentTime_(0.0),
arg_(0),
dest_(0),
spc_(spc),
gridPart_(spc_.gridPart()),
-#if WANT_DUNE_OPTIM
- linProg_( static_cast< GridViewType > (gridPart_), BoundaryValue( *this ),
- Dune::Fem::Parameter::getValue<double>("finitevolume.linearprogramming.tol", 1e-8 )
- ),
-#endif
+ linearFunction_( gridPart_ ),
+ linProg_(),
indexSet_( gridPart_.indexSet() ),
localIdSet_( gridPart_.grid().localIdSet()),
cornerPointSetContainer_(),
@@ -322,32 +415,46 @@ namespace Fem
argOrder_( spc_.order() ),
storedComboSets_(),
tolFactor_( getTolFactor() ),
- tol_1_( 1.0/getTol() ),
+ indicator_( getIndicator(parameter) ),
+ tol_1_( 2./ parameter.getValue("femdg.limiter.tolerance", double(1.0) ) ),
geoInfo_( gridPart_.indexSet() ),
faceGeoInfo_( geoInfo_.geomTypes(1) ),
phi0_( 0 ),
matrixCacheVec_( gridPart_.grid().maxLevel() + 1 ),
factors_(),
numbers_(),
- localMassMatrix_( spc_ , 2*spc_.order() ),
adaptive_((AdaptationMethodType(gridPart_.grid())).adaptive()),
cartesianGrid_( CheckCartesianType::check( gridPart_ ) ),
stepTime_(3, 0.0),
calcIndicator_(discreteModel_.calculateIndicator()),
reconstruct_(false),
- admissibleFunctions_( getAdmissibleFunctions() ),
- usedAdmissibleFunctions_( admissibleFunctions_ ),
+ admissibleFunctions_( getAdmissibleFunctions( parameter ) ),
+ usedAdmissibleFunctions_( usedAdmissibleFunctions() ),
+ extTroubledCellIndicator_( indicator_ == 3
+ ? new ModalSmoothnessIndicator< ArgumentFunctionType >() : nullptr ),
counter_( 0 )
{
- if( Parameter :: verbose () )
+ if( usedAdmissibleFunctions_ == lp )
{
- std::cout << "LimitPass: Grid is ";
+ linProg_.reset( new LinearProgramming( static_cast< GridViewType > (gridPart_), BoundaryValue( *this ),
+ parameter.getValue<double>("finitevolume.linearprogramming.tol", 1e-8 )) );
+ }
+
+ if( verbose )
+ {
+ std::cout << "LimitDGPass (Grid is ";
if( cartesianGrid_ )
- std::cout << "cartesian";
+ std::cout << "cartesian) ";
else
- std::cout << "unstructured";
- //std::cout << "! LimitEps: "<< limitEps_ << ", LimitTol: "<< 1./tol_1_ << std::endl;
- std::cout << "! Limiter.tolerance: "<< 1./tol_1_ << std::endl;
+ std::cout << "unstructured) ";
+
+ std::cout << "parameters: " << std::endl;
+ std::cout << " femdg.limiter.tolerance: " << 1./tol_1_ << std::endl;
+ std::cout << " femdg.limiter.indicator: " << indicator_ << std::endl;
+ std::string adFctName( admissibleFctNames()[ admissibleFunctions_ ] );
+ if( admissibleFunctions_ != usedAdmissibleFunctions_ )
+ adFctName += "("+admissibleFctNames()[ usedAdmissibleFunctions_ ] + ")";
+ std::cout << " femdg.limiter.admissiblefunctions: " << adFctName << std::endl;
}
// we need the flux here
@@ -356,19 +463,24 @@ namespace Fem
//! Destructor
virtual ~LimitDGPass() {
- std::cout << "~LimitDGPass: op calls " << counter_ << std::endl;
+ // std::cout << "~LimitDGPass: op calls " << counter_ << std::endl;
+ }
+
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator)
+ {
+ extTroubledCellIndicator_ = indicator;
}
//! return default face quadrature order
static int defaultVolumeQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity ));
+ return DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder( space.order( entity ) );
}
//! return default face quadrature order
static int defaultFaceQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity )) + 1;
+ return DefaultQuadrature< DiscreteFunctionSpaceType >::faceOrder( space.order( entity ) );
}
protected:
@@ -384,7 +496,6 @@ namespace Fem
return ( faceQuadOrd_ < 0 ) ? defaultFaceQuadratureOrder( spc_, entity ) : faceQuadOrd_ ;
}
-
//! get tolerance factor for shock detector
double getTolFactor() const
{
@@ -395,21 +506,43 @@ namespace Fem
}
//! get tolerance for shock detector
- double getTol() const
+ double getIndicator( const Dune::Fem::ParameterReader ¶meter ) const
{
- double tol = 1.0;
- tol = Parameter::getValue("femdg.limiter.tolerance", tol );
- return tol;
+ static const std::string indicators[] = { "user", "none", "jump" ,"modal", "always" };
+ std::string key( "femdg.limiter.indicator" );
+ return parameter.getEnum( key, indicators, 2);
+ }
+
+ static const std::string (&admissibleFctNames())[5]
+ {
+ static const std::string admissiblefct[] = { "dgonly" , "muscl", "muscl+dg", "lp", "default" };
+ return admissiblefct;
}
//! get tolerance for shock detector
- AdmissibleFunctions getAdmissibleFunctions() const
+ AdmissibleFunctions getAdmissibleFunctions( const Dune::Fem::ParameterReader ¶meter ) const
+ {
+ std::string key( "femdg.limiter.admissiblefunctions" );
+ return (AdmissibleFunctions) parameter.getEnum( key, admissibleFctNames(), _default );
+ }
+
+ AdmissibleFunctions usedAdmissibleFunctions() const
{
- // default value
- int val = 1;
- val = Parameter::getValue("femdg.limiter.admissiblefunctions", val);
- assert( val >= DGFunctions || val <= BothFunctions );
- return (AdmissibleFunctions) val;
+ if( admissibleFunctions_ == _default )
+ {
+ if( ! geoInfo_.multipleGeomTypes() && geoInfo_.geomTypes(0)[0].isSimplex() )
+ {
+ // default for simplices is muscl reconstruction because it
+ // seems to work better for simplex grids
+ return muscl;
+ }
+ else
+ {
+ // default for all other element types is lp reconstruction
+ return lp;
+ }
+ }
+ return admissibleFunctions_;
}
template <class S1, class S2>
@@ -551,7 +684,7 @@ namespace Fem
return tmp;
}
- size_t leafElements() const
+ size_t numberOfElements() const
{
return elementCounter_;
}
@@ -578,12 +711,16 @@ namespace Fem
assign( U , dest, firstThread );
// if case of finite volume scheme set admissible functions to reconstructions
- usedAdmissibleFunctions_ = reconstruct_ ? ReconstructedFunctions : admissibleFunctions_;
+ if( reconstruct_ && ! linProg_ )
+ {
+ // if linProg is not available then the only other FV reconstruction is muscl
+ usedAdmissibleFunctions_ = muscl;
+ }
// in case of reconstruction
if( reconstruct_ )
{
- // in case of non-adaptive scheme indicator not needed
+ // for FV scheme in case of non-adaptive scheme indicator not needed
calcIndicator_ = adaptive_;
// adjust quadrature orders
@@ -602,7 +739,7 @@ namespace Fem
currentTime_ = this->time();
// initialize caller
- caller_ = new DiscreteModelCallerType( *arg_, discreteModel_ );
+ caller_.reset( new DiscreteModelCallerType( *arg_, discreteModel_ ) );
caller_->setTime(currentTime_);
// calculate maximal indicator (if necessary)
@@ -624,21 +761,12 @@ namespace Fem
matrixCacheVec_.resize( numLevels );
}
-#if WANT_DUNE_OPTIM
- // helper class for evaluation of average value of discrete function
- EvalAverage average( *this, U, discreteModel_);
-
- values_.resize( size );
- // evaluate average values on all cells
- for( const auto& element : Dune::elements( gridPart_ ) )
+ if( linProg_ )
{
- average.evaluate( element, values_[ gridPart_.indexSet().index( element ) ] );
+ values_.resize( size );
+ valuesComputed_.resize( size );
+ std::fill( valuesComputed_.begin(), valuesComputed_.end(), false );
}
-
- gradients_.resize( size );
- // get reconstructions
- linProg_( gridPart_.indexSet(), values_, gradients_ );
-#endif
}
//! Some management (interface version)
@@ -667,11 +795,9 @@ namespace Fem
}
// finalize caller
- if( caller_ )
- {
- delete caller_;
- caller_ = 0;
- }
+ caller_.reset();
+ // reset usedAdmissibleFunction value to class default
+ usedAdmissibleFunctions_ = usedAdmissibleFunctions();
}
//! apply local is virtual
@@ -739,9 +865,6 @@ namespace Fem
// timer for shock detection
Dune::Timer indiTime;
- // extract types
- enum { dim = EntityType :: dimension };
-
// check argument is not zero
assert( arg_ );
@@ -749,6 +872,9 @@ namespace Fem
// set entity to caller
caller().setEntity( en );
+ // initialize linear function
+ linearFunction_.bind( en );
+
// get function to limit
const ArgumentFunctionType &U = *(std::get< argumentPosition >( *arg_ ));
@@ -756,13 +882,17 @@ namespace Fem
const LocalFunctionType uEn = U.localFunction(en);
// get geometry
- const Geometry& geo = en.geometry();
+ const Geometry& geo = linearFunction_.geometry();
// cache geometry type
const GeometryType geomType = geo.type();
// get bary center of element
- const LocalDomainType& x = geoInfo_.localCenter( geomType );
- const DomainType enBary = geo.global( x );
+ const LocalDomainType& wLocal = geoInfo_.localCenter( geomType );
+
+ DomainType& enBary = linearFunction_.center_;
+
+ // compute barycenter of element
+ enBary = geomType.isNone() ? geo.center() : geo.global( wLocal );
const IntersectionIteratorType endnit = gridPart_.iend( en );
IntersectionIteratorType niter = gridPart_.ibegin(en);
@@ -774,7 +904,7 @@ namespace Fem
RangeType shockIndicator(0);
RangeType adaptIndicator(0);
- RangeType enVal;
+ RangeType& enVal = linearFunction_.value_;
// if limiter is true then limitation is done
// when we want to reconstruct in any case then
@@ -796,9 +926,8 @@ namespace Fem
// otherwise everything may be corrupted
{
// get barycenter of entity
- const DomainType& enBaryLocal = (int(dimGrid) == int(dimDomain)) ?
- geoInfo_.localCenter( geomType ) :
- geo.local( enBary ) ;
+ const DomainType& enBaryLocal = (int(dimGrid) == int(dimDomain) && ! geomType.isNone() ) ?
+ wLocal : geo.local( enBary ) ;
// check that average value is physical
if( discreteModel_.hasPhysical() &&
@@ -824,8 +953,8 @@ namespace Fem
else if ( calcIndicator_ ) // otherwise compute shock indicator
{
// check shock indicator
- limiter = calculateIndicator(en, uEn, geo, limiter, // parameter
- limit, shockIndicator, adaptIndicator); // return values
+ limiter = calculateIndicator(U, uEn, geo, limiter, // parameter
+ limit, shockIndicator, adaptIndicator); // return values
}
else if( !reconstruct_ )
{
@@ -879,10 +1008,12 @@ namespace Fem
EvalAverage average( *this, U, discreteModel_, geo.volume() );
// setup neighbors barycenter and mean value for all neighbors
- LimiterUtilityType::setupNeighborValues( gridPart_, en, average, enBary, enVal,
+ LimiterUtilityType::setupNeighborValues( gridPart_, en, average, enBary, enVal, uEn,
StructuredGrid, flags, barys_, nbVals_ );
}
+ const unsigned int enIndex = indexSet_.index( en );
+
// mark entity as finished, even if not limited everything necessary was done
assert( indexSet_.index( en ) < int(visited_.size()) );
visited_[ indexSet_.index( en ) ] = true ;
@@ -893,15 +1024,18 @@ namespace Fem
// increase number of limited elements
++limitedElements_;
- const unsigned int enIndex = indexSet_.index( en );
-#if WANT_DUNE_OPTIM
- bool useLinProg = true ;
- if( useLinProg )
+ GradientType& gradient = linearFunction_.gradient_;
+
+ // use linear function from LP reconstruction
+ if( usedAdmissibleFunctions_ == lp )
{
- deoMod_ = gradients_[ enIndex ];
+ // compute average values on neighboring elements
+ fillAverageValues( en, enIndex, U, enVal );
+ assert( linProg_ );
+ // compute optimal linear reconstruction
+ linProg_->applyLocal( en, gridPart_.indexSet(), values_, gradient );
}
else
-#endif
{
// obtain combination set
ComboSetType& comboSet = storedComboSets_[ nbVals_.size() ];
@@ -912,26 +1046,26 @@ namespace Fem
}
// reset values
- deoMods_.clear();
+ gradients_.clear();
comboVec_.clear();
- if( usedAdmissibleFunctions_ >= ReconstructedFunctions )
+ if( usedAdmissibleFunctions_ == muscl || usedAdmissibleFunctions_ == muscldg )
{
// level is only needed for Cartesian grids to access the matrix caches
const int matrixCacheLevel = ( flags.cartesian ) ? en.level() : 0 ;
assert( matrixCacheLevel < (int) matrixCacheVec_.size() );
MatrixCacheType& matrixCache = matrixCacheVec_[ matrixCacheLevel ];
- // calculate linear functions, stored in deoMods_ and comboVec_
+ // calculate linear functions, stored in gradients_ and comboVec_
LimiterUtilityType::calculateLinearFunctions( comboSet, geomType, flags,
barys_, nbVals_,
matrixCache,
- deoMods_,
+ gradients_,
comboVec_ );
}
// add DG Function
- if( (usedAdmissibleFunctions_ % 2) == DGFunctions )
+ if( usedAdmissibleFunctions_ == dgonly || usedAdmissibleFunctions_ == muscldg )
{
addDGFunction( en, geo, uEn, enVal, enBary );
}
@@ -939,10 +1073,10 @@ namespace Fem
// Limiting
std::vector< RangeType > factors;
LimiterUtilityType::limitFunctions(
- discreteModel_.limiterFunction(), comboVec_, barys_, nbVals_, deoMods_, factors );
+ discreteModel_.limiterFunction(), comboVec_, barys_, nbVals_, gradients_, factors );
// take maximum of limited functions
- LimiterUtilityType::getMaxFunction(deoMods_, deoMod_, factors_[ enIndex ], numbers_[ enIndex ], factors );
+ LimiterUtilityType::getMaxFunction(gradients_, gradient, factors_[ enIndex ], numbers_[ enIndex ], factors );
} // end if linProg
////////////////////////////////////////////////////////////////////
@@ -954,8 +1088,8 @@ namespace Fem
// get local funnction for limited values
DestLocalFunctionType limitEn = dest_->localFunction(en);
- // project deoMod_ to limitEn
- L2project(en, geo, enBary, enVal, limit, deoMod_, limitEn );
+ // project linearFunction onto limitEn
+ interpolate( en, geo, limit, linearFunction_, limitEn );
assert( checkPhysical(en, geo, limitEn) );
@@ -965,6 +1099,34 @@ namespace Fem
}
protected:
+ void fillAverageValues( const EntityType& en, const unsigned int enIndex,
+ const ArgumentFunctionType &U, const RangeType& enVal ) const
+ {
+ // helper class for evaluation of average value of discrete function
+ EvalAverage average( *this, U, discreteModel_);
+
+ if( ! valuesComputed_[ enIndex ] )
+ {
+ values_[ enIndex ] = enVal;
+ valuesComputed_[ enIndex ] = true;
+ }
+
+ const auto& indexSet = gridPart_.indexSet();
+ for (const auto& intersection : intersections(gridPart_, en) )
+ {
+ if( intersection.neighbor() )
+ {
+ const auto neighbor = intersection.outside();
+ const unsigned int nbIndex = indexSet.index( neighbor );
+ if( ! valuesComputed_[ nbIndex ] )
+ {
+ average.evaluate( neighbor, values_[ nbIndex ] );
+ valuesComputed_[ nbIndex ] = true;
+ }
+ }
+ }
+ }
+
// add linear components of the DG function
void addDGFunction(const EntityType& en,
const Geometry& geo,
@@ -1014,7 +1176,7 @@ namespace Fem
A.mv( rhs, D[ r ]);
}
- deoMods_.push_back( D );
+ gradients_.push_back( D );
std::vector<int> comb( nbVals_.size() );
const size_t combSize = comb.size();
for (size_t i=0;i<combSize; ++i) comb[ i ] = i;
@@ -1050,7 +1212,6 @@ namespace Fem
const Geometry& geo,
const LocalFunctionImp& uEn) const
{
- enum { dim = dimGrid };
if( discreteModel_.hasPhysical() )
{
if( en.type().isNone() )
@@ -1104,63 +1265,42 @@ namespace Fem
}
};
+
// L2 projection
template <class LocalFunctionImp>
- void L2project(const EntityType& en,
- const Geometry& geo,
- const DomainType& enBary,
- const RangeType& enVal,
- const FieldVector<bool,dimRange>& limit,
- const GradientType& deoMod,
- LocalFunctionImp& limitEn,
- const bool constantValue = false ) const
+ void interpolate(const EntityType& en,
+ const Geometry& geo,
+ const FieldVector<bool,dimRange>& limit,
+ const LinearFunction& linearFunction,
+ LocalFunctionImp& limitEn ) const
{
- enum { dim = dimGrid };
-
// true if geometry mapping is affine
- const bool affineMapping = localMassMatrix_.affine();
+ const bool affineMapping = geo.affine();
// set zero dof to zero
uTmpLocal_.init( en );
uTmpLocal_.clear();
- // get quadrature for destination space order
- VolumeQuadratureType quad( en, spc_.order() + 1 );
+ const auto interpolation = spc_.interpolation( en );
- const int quadNop = quad.nop();
- for(int qP = 0; qP < quadNop ; ++qP)
+ const bool constantValue = linearFunction.isConstant();
+ if( constantValue )
{
- // get quadrature weight
- const double intel = (affineMapping) ?
- quad.weight(qP) : // affine case
- quad.weight(qP) * geo.integrationElement( quad.point(qP) ); // general case
-
- RangeType retVal( enVal );
- // if we don't have only constant value then evaluate function
- if( !constantValue )
- {
- // get global coordinate
- DomainType point = geo.global( quad.point( qP ) );
- point -= enBary;
-
- // evaluate linear function
- for( int r = 0; r < dimRange; ++r )
- retVal[ r ] += (deoMod[ r ] * point);
- }
-
- retVal *= intel;
- uTmpLocal_.axpy( quad[ qP ], retVal );
+ const ConstantFunction& constFct = linearFunction;
+ interpolation( constFct, uTmpLocal_.localDofVector() );
+ }
+ else
+ {
+ interpolation( linearFunction, uTmpLocal_.localDofVector() );
}
-
- if( !affineMapping )
- localMassMatrix_.applyInverse( en, uTmpLocal_ );
// check physicality of projected data
if ( (! constantValue) && (! checkPhysical(en, geo, uTmpLocal_)) )
{
// for affine mapping we only need to set higher moments to zero
- if( affineMapping )
- {
+ if( Dune::Fem::Capabilities::isHierarchic< DiscreteFunctionSpaceType > :: v &&
+ affineMapping )
+ { // note: the following assumes an ONB made up of moments and affine mapping
const int numBasis = uTmpLocal_.numDofs()/dimRange;
for(int i=1; i<numBasis; ++i)
{
@@ -1173,9 +1313,9 @@ namespace Fem
}
else
{
- // if not physical project to mean value
- L2project(en,geo,enBary,enVal,limit,deoMod_, limitEn, true);
- return ;
+ // general interpolation of constant value
+ const ConstantFunction& constFct = linearFunction;
+ interpolation( constFct, uTmpLocal_.localDofVector() );
}
}
@@ -1189,7 +1329,7 @@ namespace Fem
{
// copy to limitEn skipping components that should not be limited
const int numBasis = uTmpLocal_.numDofs()/dimRange;
- for(int i=1; i<numBasis; ++i)
+ for(int i=0; i<numBasis; ++i)
{
for( const auto& r : discreteModel_.model().limitedRange() )
{
@@ -1202,7 +1342,7 @@ namespace Fem
template <class BasisFunctionSetType, class PointType>
const RangeType& evaluateConstantBasis( const BasisFunctionSetType& basisSet,
- const PointType& x ) const
+ const PointType& x ) const
{
// calculate constant part of the basis functions
if( ! (phi0_[ 0 ] > 0 ) )
@@ -1231,8 +1371,9 @@ namespace Fem
RangeType& val) const
{
bool notphysical = false;
- if( Dune::Fem::Capabilities::isHierarchic< DiscreteFunctionSpaceType > :: v
- && localMassMatrix_.affine() )
+ const Geometry& geo = en.geometry();
+
+ if( Dune::Fem::Capabilities::isHierarchic< DiscreteFunctionSpaceType > :: v && geo.affine() )
{
// get point quadrature
VolumeQuadratureType quad( en, 0 );
@@ -1253,8 +1394,6 @@ namespace Fem
}
else
{
- const Geometry& geo = en.geometry();
-
// get quadrature
VolumeQuadratureType quad( en, volumeQuadratureOrder( en ) );
@@ -1287,11 +1426,22 @@ namespace Fem
return notphysical;
}
- template <bool conforming>
bool integrateIntersection(const IntersectionType & intersection,
const EntityType& nb,
RangeType& shockIndicator,
RangeType& adaptIndicator) const
+ {
+ if( ! conformingGridPart && ! intersection.conforming() )
+ return integrateIntersectionImpl< false >(intersection, nb, shockIndicator, adaptIndicator);
+ else
+ return integrateIntersectionImpl< true >(intersection, nb, shockIndicator, adaptIndicator);
+ }
+
+ template <bool conforming>
+ bool integrateIntersectionImpl(const IntersectionType & intersection,
+ const EntityType& nb,
+ RangeType& shockIndicator,
+ RangeType& adaptIndicator) const
{
// make sure we got the right conforming statement
assert( intersection.conforming() == conforming );
@@ -1381,15 +1531,15 @@ namespace Fem
}
// calculate shock detector
- bool calculateIndicator(const EntityType& en,
+ bool calculateIndicator(const ArgumentFunctionType& U,
const LocalFunctionType& uEn,
const Geometry& geo,
const bool initLimiter,
FieldVector<bool,dimRange>& limit,
RangeType& shockIndicator,
- RangeType& adaptIndicator) const
+ RangeType& adaptIndicator ) const
{
- enum { dim = EntityType :: dimension };
+ const EntityType& en = uEn.entity();
// calculate circume during neighbor check
double circume = 0.0;
@@ -1442,7 +1592,6 @@ namespace Fem
shockIndicator = -1;
return true;
}
-
}
else
{ // non-conforming case
@@ -1496,52 +1645,46 @@ namespace Fem
// add face vol to circume
circume += vol;
- // order of quadrature
- const int jumpQuadOrd = spc_.order();
-
- // check all neighbors
- if (intersection.neighbor())
+ // internal shock indicator based on Krivodonova et al.
+ if( ! extTroubledCellIndicator_ )
{
- // get neighbor entity
- const EntityType& nb = intersection.outside();
+ // order of quadrature
+ const int jumpQuadOrd = spc_.order();
- // conforming case
- if( ! conformingGridPart && ! intersection.conforming() )
- { // non-conforming case
- if (integrateIntersection< false > (intersection, nb, shockIndicator, adaptIndicator))
+ // check all neighbors
+ if (intersection.neighbor())
+ {
+ // get neighbor entity
+ const EntityType& nb = intersection.outside();
+
+ if (integrateIntersection(intersection, nb, shockIndicator, adaptIndicator))
{
shockIndicator = -1;
return true;
}
}
- else
+
+ // check all neighbors
+ if ( intersection.boundary() )
{
- if (integrateIntersection< true > (intersection, nb, shockIndicator, adaptIndicator))
+ FaceQuadratureType faceQuadInner(gridPart_,intersection, jumpQuadOrd, FaceQuadratureType::INSIDE);
+
+ // initialize intersection
+ caller().initializeBoundary( intersection, faceQuadInner );
+
+ if (integrateBoundary(en, intersection, faceQuadInner, shockIndicator, adaptIndicator))
{
shockIndicator = -1;
return true;
}
}
}
-
- // check all neighbors
- if ( intersection.boundary() )
- {
- FaceQuadratureType faceQuadInner(gridPart_,intersection, jumpQuadOrd, FaceQuadratureType::INSIDE);
-
- // initialize intersection
- caller().initializeBoundary( intersection, faceQuadInner );
-
- if (integrateBoundary(en, intersection, faceQuadInner, shockIndicator, adaptIndicator))
- {
- shockIndicator = -1;
- return true;
- }
- }
}
-
} // end intersection iterator
+ if (indicator_ == 1)
+ return false;
+
// calculate max face volume
{
// min faceVol
@@ -1566,11 +1709,34 @@ namespace Fem
// multiply h pol ord with circume
const double circFactor = (circume > 0.0) ? (hPowPolOrder/(circume * tolFactor_ )) : 0.0;
+ assert(indicator_!=0 || extTroubledCellIndicator_);
+ if( extTroubledCellIndicator_ ) // also true for indicator_=3
+ {
+ shockIndicator = ( tol_1_ * (*extTroubledCellIndicator_)( U, uEn ) );
+ }
+ else if( indicator_ == 2 )
+ {
+ for(int r=0; r<dimRange; ++r)
+ {
+ shockIndicator[r] = std::abs(shockIndicator[r]) * circFactor * tol_1_;
+ }
+ }
+ else if ( indicator_ == 4 )
+ {
+ // always limit on each cell
+ shockIndicator = 2.;
+ }
+ else
+ {
+ std::cout << "wrong indicator selection\n";
+ assert(0);
+ abort();
+ }
+
for(int r=0; r<dimRange; ++r)
{
- // only scale shock indicator with tolerance
- shockIndicator[r] = std::abs(shockIndicator[r]) * circFactor * tol_1_;
adaptIndicator[r] = std::abs(adaptIndicator[r]) * circFactor;
+
if(shockIndicator[r] > 1.)
{
limit[r] = true;
@@ -1589,10 +1755,10 @@ namespace Fem
{
// check whether we have an inflow intersection or not
const int quadNop = quad.nop();
- for(int l=0; l<quadNop; ++l)
+ for(int qp=0; qp<quadNop; ++qp)
{
// check physicality of value
- const bool physical = caller().checkPhysical( intersection, quad, l );
+ const bool physical = caller().checkPhysical( intersection, quad, qp );
if( checkPhysical && ! physical )
{
@@ -1603,7 +1769,7 @@ namespace Fem
if( ! inflowIntersection )
{
// check intersection
- if( physical && caller().checkDirection(intersection, quad, l) )
+ if( physical && caller().checkDirection(intersection, quad, qp) )
{
inflowIntersection = true;
// in case of physicality check is disabled
@@ -1634,7 +1800,7 @@ namespace Fem
}
private:
- mutable DiscreteModelCallerType *caller_;
+ mutable std::unique_ptr< DiscreteModelCallerType > caller_;
DiscreteModelType& discreteModel_;
mutable double currentTime_;
@@ -1644,9 +1810,9 @@ namespace Fem
const DiscreteFunctionSpaceType& spc_;
GridPartType& gridPart_;
-#if WANT_DUNE_OPTIM
- mutable LinearProgramming linProg_;
-#endif
+ mutable LinearFunction linearFunction_;
+
+ std::unique_ptr< LinearProgramming > linProg_;
const IndexSetType& indexSet_;
const LocalIdSetType& localIdSet_;
@@ -1665,16 +1831,16 @@ namespace Fem
// tolerance to scale shock indicator
const double tolFactor_;
+ std::size_t indicator_;
const double tol_1_;
// if true scheme is TVD
const GeometryInformationType geoInfo_;
const FaceGeometryInformationType faceGeoInfo_;
- mutable GradientType deoMod_;
mutable RangeType phi0_ ;
- mutable std::vector< GradientType > deoMods_;
+ mutable std::vector< GradientType > gradients_;
mutable std::vector< CheckType > comboVec_;
mutable std::vector< DomainType > barys_;
@@ -1687,7 +1853,6 @@ namespace Fem
// vector for stroing the information which elements have been computed already
mutable std::vector< bool > visited_;
- LocalMassMatrixType localMassMatrix_;
//! true if limiter is used in adaptive scheme
const bool adaptive_;
//! true if grid is cartesian like
@@ -1696,18 +1861,21 @@ namespace Fem
mutable std::vector<double> stepTime_;
mutable size_t elementCounter_;
- //! true if indicator should be calculated
+ //! true if troubled cell indicator should be calculated
mutable bool calcIndicator_;
//! true if limiter is used as finite volume scheme of higher order
mutable bool reconstruct_;
// choice of admissible linear functions
- const AdmissibleFunctions admissibleFunctions_;
- mutable AdmissibleFunctions usedAdmissibleFunctions_ ;
+ const AdmissibleFunctions admissibleFunctions_;
+ mutable AdmissibleFunctions usedAdmissibleFunctions_ ;
mutable std::vector< RangeType > values_;
- mutable std::vector< GradientType > gradients_;
+ mutable std::vector< bool > valuesComputed_;
+
+ // function pointer to external troubled cell indicator, can be nullptr
+ TroubledCellIndicatorType extTroubledCellIndicator_;
mutable int counter_;
}; // end DGLimitPass
diff --git a/dune/fem-dg/operator/limiter/lpreconstruction.hh b/dune/fem-dg/operator/limiter/lpreconstruction.hh
index fff7baeac5fb6882b25459163fca823da9098ebd..c6269e30ce9d81d9a1993506d43e00234b5f70b7 100644
--- a/dune/fem-dg/operator/limiter/lpreconstruction.hh
+++ b/dune/fem-dg/operator/limiter/lpreconstruction.hh
@@ -83,6 +83,7 @@ namespace Dune
typedef Optim::GaussJordanSolver< FieldMatrix< Field, GlobalCoordinate::dimension, GlobalCoordinate::dimension > > Solver;
typedef Optim::LinearProgramming< Solver, false > LP;
+ typedef std::vector< std::pair< GlobalCoordinate, StateVector > > DifferencesVectorType;
public:
LPReconstruction ( const GridView &gridView, BoundaryValue boundaryValue, Real tolerance )
: gridView_( gridView ),
@@ -105,119 +106,127 @@ namespace Dune
}
}
- template< class Mapper, class Vector >
- void operator () ( const Mapper &mapper, const Vector &u, std::vector< Jacobian > &du ) const
+ template< class Entity, class Mapper, class Vector >
+ void applyLocal ( const Entity& element,
+ const Mapper &mapper,
+ const Vector &u,
+ Jacobian& du ) const
{
- du.resize( u.size() );
+ DifferencesVectorType& differences = differences_;
+ differences.clear();
- std::vector< std::pair< GlobalCoordinate, StateVector > > differences;
Constraints constraints;
- const auto end = gridView().template end< 0, Dune::InteriorBorder_Partition >();
- for( auto it = gridView().template begin< 0, Dune::InteriorBorder_Partition>(); it != end; ++it )
+ const std::size_t elIndex = mapper.index( element );
+ const GlobalCoordinate elCenter = element.geometry().center();
+
+ std::array< unsigned int, dimension+1 > select;
+ const std::vector< unsigned int > &faceAxes = faceAxes_[ LocalGeometryTypeIndex::index( element.type() ) ];
+ if( !faceAxes.empty() )
{
- const auto element = *it;
+ const auto iend = gridView().iend( element );
+ for( auto iit = gridView().ibegin( element ); iit != iend; ++iit )
+ {
+ const auto intersection = *iit;
- const std::size_t elIndex = mapper.index( element );
- const GlobalCoordinate elCenter = element.geometry().center();
+ select[ faceAxes[ intersection.indexInInside() ] ] = differences.size();
- std::array< unsigned int, dimension+1 > select;
- const std::vector< unsigned int > &faceAxes = faceAxes_[ LocalGeometryTypeIndex::index( element.type() ) ];
- if( !faceAxes.empty() )
- {
- differences.clear();
- const auto iend = gridView().iend( element );
- for( auto iit = gridView().ibegin( element ); iit != iend; ++iit )
+ if( intersection.boundary() )
+ {
+ const GlobalCoordinate iCenter = intersection.geometry().center();
+ const GlobalCoordinate iNormal = intersection.centerUnitOuterNormal();
+ const StateVector uBnd = boundaryValue_( intersection, iCenter, iNormal, u[ elIndex ] );
+ differences.emplace_back( iCenter - elCenter, uBnd - u[ elIndex ] );
+ }
+ else if( intersection.neighbor() )
{
- const auto intersection = *iit;
-
- select[ faceAxes[ intersection.indexInInside() ] ] = differences.size();
-
- if( intersection.boundary() )
- {
- const GlobalCoordinate iCenter = intersection.geometry().center();
- const GlobalCoordinate iNormal = intersection.centerUnitOuterNormal();
- const StateVector uBnd = boundaryValue_( intersection, iCenter, iNormal, u[ elIndex ] );
- differences.emplace_back( iCenter - elCenter, uBnd - u[ elIndex ] );
- }
- else if( intersection.neighbor() )
- {
- const auto neighbor = intersection.outside();
- const GlobalCoordinate nbCenter = neighbor.geometry().center();
- differences.emplace_back( nbCenter - elCenter, u[ mapper.index( neighbor ) ] - u[ elIndex ] );
- }
+ const auto neighbor = intersection.outside();
+ const GlobalCoordinate nbCenter = neighbor.geometry().center();
+ differences.emplace_back( nbCenter - elCenter, u[ mapper.index( neighbor ) ] - u[ elIndex ] );
}
}
- else
+ }
+ else
+ {
+ Dune::ReservedVector< GlobalCoordinate, dimension > onb;
+
+ const auto iend = gridView().iend( element );
+ for( auto iit = gridView().ibegin( element ); iit != iend; ++iit )
{
- Dune::ReservedVector< GlobalCoordinate, dimension > onb;
+ const auto intersection = *iit;
+
+ select[ onb.size() ] = differences.size();
- differences.clear();
- const auto iend = gridView().iend( element );
- for( auto iit = gridView().ibegin( element ); iit != iend; ++iit )
+ if( intersection.boundary() )
{
- const auto intersection = *iit;
-
- select[ onb.size() ] = differences.size();
-
- if( intersection.boundary() )
- {
- const GlobalCoordinate iCenter = intersection.geometry().center();
- const GlobalCoordinate iNormal = intersection.centerUnitOuterNormal();
- const StateVector uBnd = boundaryValue_( intersection, iCenter, iNormal, u[ elIndex ] );
- differences.emplace_back( iCenter - elCenter, uBnd - u[ elIndex ] );
- }
- else if( intersection.neighbor() )
- {
- const auto neighbor = intersection.outside();
- const GlobalCoordinate nbCenter = neighbor.geometry().center();
- differences.emplace_back( nbCenter - elCenter, u[ mapper.index( neighbor ) ] - u[ elIndex ] );
- }
-
- if( onb.size() < dimension )
- {
- GlobalCoordinate dx = differences.back().first;
- for( const GlobalCoordinate &v : onb )
- dx.axpy( -(dx*v), v );
-
- const Real dxNorm = dx.two_norm();
- if( dxNorm >= tolerance_ )
- onb.push_back( dx /= dxNorm );
- }
+ const GlobalCoordinate iCenter = intersection.geometry().center();
+ const GlobalCoordinate iNormal = intersection.centerUnitOuterNormal();
+ const StateVector uBnd = boundaryValue_( intersection, iCenter, iNormal, u[ elIndex ] );
+ differences.emplace_back( iCenter - elCenter, uBnd - u[ elIndex ] );
+ }
+ else if( intersection.neighbor() )
+ {
+ const auto neighbor = intersection.outside();
+ const GlobalCoordinate nbCenter = neighbor.geometry().center();
+ differences.emplace_back( nbCenter - elCenter, u[ mapper.index( neighbor ) ] - u[ elIndex ] );
}
- }
- // reserve memory for constraints
- const std::size_t numConstraints = differences.size();
- constraints.resize( 2u*numConstraints );
- Optim::ActiveIndexMapper< SmallObjectAllocator< unsigned int > > active( GlobalCoordinate::dimension, constraints.size() );
- for( int j = 0; j < StateVector::dimension; ++j )
- {
- GlobalCoordinate negGradient( 0 );
- for( std::size_t i = 0u; i < numConstraints; ++i )
+ if( onb.size() < dimension )
{
- const Field sign = (differences[ i ].second[ j ] >= 0 ? 1 : -1);
+ GlobalCoordinate dx = differences.back().first;
+ for( const GlobalCoordinate &v : onb )
+ dx.axpy( -(dx*v), v );
- negGradient.axpy( sign, differences[ i ].first );
+ const Real dxNorm = dx.two_norm();
+ if( dxNorm >= tolerance_ )
+ onb.push_back( dx /= dxNorm );
+ }
+ }
+ }
- constraints[ 2*i ].normal() = differences[ i ].first;
- constraints[ 2*i ].normal() *= sign;
- constraints[ 2*i ].rhs() = sign*differences[ i ].second[ j ];
+ // reserve memory for constraints
+ const std::size_t numConstraints = differences.size();
+ constraints.resize( 2u*numConstraints );
+ Optim::ActiveIndexMapper< SmallObjectAllocator< unsigned int > > active( GlobalCoordinate::dimension, constraints.size() );
+ for( int j = 0; j < StateVector::dimension; ++j )
+ {
+ GlobalCoordinate negGradient( 0 );
+ for( std::size_t i = 0u; i < numConstraints; ++i )
+ {
+ const Field sign = (differences[ i ].second[ j ] >= 0 ? 1 : -1);
- constraints[ 2*i+1 ].normal() = constraints[ 2*i ].normal();
- constraints[ 2*i+1 ].normal() *= Field( -1 );
- constraints[ 2*i+1 ].rhs() = 0;
- }
+ negGradient.axpy( sign, differences[ i ].first );
- // activate GlobalCoordinate::dimension constraints active in the origin
- active.clear();
- for( int i = 0; i < dimension; ++i )
- active.activate( 2*select[ i ]+1 );
+ constraints[ 2*i ].normal() = differences[ i ].first;
+ constraints[ 2*i ].normal() *= sign;
+ constraints[ 2*i ].rhs() = sign*differences[ i ].second[ j ];
- // solve
- du[ elIndex ][ j ] = 0;
- lp_( negGradient, constraints, du[ elIndex ][ j ], active );
+ constraints[ 2*i+1 ].normal() = constraints[ 2*i ].normal();
+ constraints[ 2*i+1 ].normal() *= Field( -1 );
+ constraints[ 2*i+1 ].rhs() = 0;
}
+
+ // activate GlobalCoordinate::dimension constraints active in the origin
+ active.clear();
+ for( int i = 0; i < dimension; ++i )
+ active.activate( 2*select[ i ]+1 );
+
+ // solve
+ du[ j ] = 0;
+ lp_( negGradient, constraints, du[ j ], active );
+ }
+ }
+
+ template< class Mapper, class Vector >
+ void operator () ( const Mapper &mapper, const Vector &u, std::vector< Jacobian > &du ) const
+ {
+ du.resize( u.size() );
+
+ const auto end = gridView().template end< 0, Dune::InteriorBorder_Partition >();
+ for( auto it = gridView().template begin< 0, Dune::InteriorBorder_Partition>(); it != end; ++it )
+ {
+ const auto element = *it;
+ applyLocal( element, mapper, u, du[ mapper.index( element ) ] );
}
//auto handle = vectorCommDataHandle( mapper, du, [] ( Jacobian a, Jacobian b ) { return b; } );
@@ -232,6 +241,7 @@ namespace Dune
Real tolerance_;
LP lp_;
std::vector< std::vector< unsigned int > > faceAxes_;
+ mutable DifferencesVectorType differences_;
};
diff --git a/dune/fem-dg/operator/limiter/scalinglimitpass.hh b/dune/fem-dg/operator/limiter/scalinglimitpass.hh
index 9f0825ac4510f4a6cbf9af406bb7aba77f877153..c4d1d8d7934d78c79d388f934c0ae208b5405f13 100644
--- a/dune/fem-dg/operator/limiter/scalinglimitpass.hh
+++ b/dune/fem-dg/operator/limiter/scalinglimitpass.hh
@@ -116,24 +116,9 @@ namespace Fem
//! type of cartesian grid checker
typedef CheckCartesian< GridPartType > CheckCartesianType;
- protected:
- template <class DiscreteSpace>
- struct HierarchicalBasis
- {
- static const bool v = false ;
- };
-
- template < class FunctionSpace, class GridPart, int polOrder, template< class > class Storage >
- struct HierarchicalBasis< DiscontinuousGalerkinSpace< FunctionSpace, GridPart, polOrder, Storage > >
- {
- static const bool v = true ;
- };
+ //! function describing an external troubled cell indicator
+ typedef std::shared_ptr< TroubledCellIndicatorBase<ArgumentFunctionType> > TroubledCellIndicatorType;
- template < class FunctionSpace, class GridPart, int polOrder, template< class > class Storage >
- struct HierarchicalBasis< HierarchicLegendreDiscontinuousGalerkinSpace< FunctionSpace, GridPart, polOrder, Storage > >
- {
- static const bool v = true ;
- };
public:
//- Public methods
@@ -146,10 +131,30 @@ namespace Fem
* \param fQ order of face quadrature
*/
ScalingLimitDGPass(DiscreteModelType& problem,
- PreviousPassType& pass,
- const DiscreteFunctionSpaceType& spc,
- const int vQ = -1,
- const int fQ = -1 ) :
+ PreviousPassType& pass,
+ const DiscreteFunctionSpaceType& spc,
+ const int vQ = -1,
+ const int fQ = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose()) :
+ ScalingLimitDGPass( problem, pass, spc, Dune::Fem::Parameter::container() )
+ {}
+
+ //- Public methods
+ /** \brief constructor
+ *
+ * \param problem Actual problem definition (see problem.hh)
+ * \param pass Previous pass
+ * \param spc Space belonging to the discrete function local to this pass
+ * \param vQ order of volume quadrature
+ * \param fQ order of face quadrature
+ */
+ ScalingLimitDGPass(DiscreteModelType& problem,
+ PreviousPassType& pass,
+ const DiscreteFunctionSpaceType& spc,
+ const Dune::Fem::ParameterReader ¶meter = Dune::Fem::Parameter::container(),
+ const int vQ = -1,
+ const int fQ = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose()) :
BaseType(pass, spc),
caller_( 0 ),
discreteModel_(problem),
@@ -158,6 +163,7 @@ namespace Fem
dest_(0),
spc_(spc),
gridPart_(spc_.gridPart()),
+ scaledFunction_( gridPart_, discreteModel_.model(), spc_.order() ),
indexSet_( gridPart_.indexSet() ),
cornerPointSetContainer_(),
dofConversion_(dimRange),
@@ -166,7 +172,7 @@ namespace Fem
argOrder_( spc_.order() ),
geoInfo_( gridPart_.indexSet() ),
phi0_( 0 ),
- localMassMatrix_( spc_ , 2*spc_.order()+3 ),
+ localMassMatrix_( spc_, [this](const int order) { return DefaultQuadrature<DiscreteFunctionSpaceType >::volumeOrder(order); } ),
stepTime_(3, 0.0)
{
// we need the flux here
@@ -182,31 +188,25 @@ namespace Fem
//! return default face quadrature order
static int defaultVolumeQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity ));
+ return DefaultQuadrature<DiscreteFunctionSpaceType >::volumeOrder( space.order( entity ) );
}
//! return default face quadrature order
static int defaultFaceQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity )) + 1;
+ return DefaultQuadrature<DiscreteFunctionSpaceType >::faceOrder( space.order( entity ) );
}
//! return appropriate quadrature order, default is 2 * order(entity)
int volumeQuadratureOrder( const EntityType& entity ) const
{
- return ( volumeQuadOrd_ < 0 ) ? ( 2*spc_.order( entity )+3 ) : volumeQuadOrd_ ;
- }
-
- //! return default face quadrature order
- int defaultFaceQuadOrder( const EntityType& entity ) const
- {
- return (2 * spc_.order( entity )) + 3;
+ return ( volumeQuadOrd_ < 0 ) ? defaultVolumeQuadratureOrder( spc_, entity ) : volumeQuadOrd_ ;
}
//! return appropriate quadrature order, default is 2 * order( entity ) + 1
int faceQuadratureOrder( const EntityType& entity ) const
{
- return ( faceQuadOrd_ < 0 ) ? defaultFaceQuadOrder( entity ) : faceQuadOrd_ ;
+ return ( faceQuadOrd_ < 0 ) ? defaultFaceQuadratureOrder( spc_, entity ) : faceQuadOrd_ ;
}
template <class S1, class S2>
@@ -339,6 +339,117 @@ namespace Fem
}
};
+
+ struct ScaledFunction :
+ public Dune::Fem::BindableGridFunction< GridPartType, Dune::Dim<dimRange> >
+ {
+ typedef Dune::Fem::BindableGridFunction<GridPartType, Dune::Dim<dimRange> > Base;
+ typedef typename Base::RangeType RangeType;
+ typedef typename Base::EntityType EntityType;
+ using Base::Base;
+
+ typedef typename DiscreteModelType :: ModelType ModelType;
+
+ protected:
+ const ModelType& model_;
+ int order_;
+
+ mutable unsigned int quadId_;
+ mutable std::vector< RangeType > tmpVal_ ;
+ mutable RangeType theta_;
+ public:
+ using Base::bind;
+ mutable RangeType enVal_;
+ mutable int count_ ;
+
+ ScaledFunction( const GridPartType& gridPart,
+ const ModelType& model, const int order )
+ : Base( gridPart ),
+ model_( model ),
+ order_( order ),
+ quadId_( -1 ),
+ tmpVal_(), theta_(1), enVal_(0), count_(0)
+ {}
+
+ void bind( const EntityType& entity, const int order )
+ {
+ Base::bind( entity );
+ order_ = order;
+ }
+
+ template <class Quadrature>
+ std::vector< RangeType >& getValues( const Quadrature& quadrature ) const
+ {
+ // adjust size if necessary
+ tmpVal_.resize( quadrature.nop() );
+
+ // store quadrature id for later check during evaluate
+ quadId_ = quadrature.id();
+ count_ = 0;
+ return tmpVal_;
+ }
+
+ // return reference to scaling factor
+ RangeType& resetTheta () const
+ {
+ // reset theta to default value which does no scaling
+ theta_ = 1;
+ return theta_;
+ }
+
+ template <class Quadrature, class RangeArray>
+ void evaluateQuadrature(const Quadrature& quadrature, RangeArray &values) const
+ {
+ assert( quadrature.id() == quadId_ );
+ const int quadNop = quadrature.nop();
+ assert( quadNop == int(tmpVal_.size()) );
+
+ for(int qP = 0; qP < quadNop ; ++qP)
+ {
+ evaluate( qP, values[ qP ] );
+ }
+ }
+
+ template <class Point>
+ void evaluate(const Point &x, typename Base::RangeType &value) const
+ {
+ evaluate( count_, value );
+ ++count_;
+ return ;
+
+ std::cout << "Wrong quadrature " << std::endl;
+ assert( false );
+ std::abort();
+ }
+
+ template <class Quadrature>
+ void evaluate(const QuadraturePointWrapper< Quadrature > &x, typename Base::RangeType &value) const
+ {
+ assert( x.quadrature().id() == quadId_ );
+ evaluate( x.index(), value );
+ }
+
+ unsigned int order() const { return order_; }
+ std::string name() const { return "ScalingLimmitPass::ScaledFunction"; } // needed for output
+
+ protected:
+ void evaluate(const int qP, typename Base::RangeType &value) const
+ {
+ assert( qP < int(tmpVal_.size() ));
+ // copy value
+ value = tmpVal_[ qP ];
+
+ // \tilde{p}(x) = \theta (p(x) - \bar{u}) + \bar{u}
+ // limitedRange contains all components that should be modified
+ // default is 0,...,dimRange-1
+ for( const auto& d : model_.limitedRange() )
+ {
+ value[ d ] = theta_[ d ] * ( value[ d ] - enVal_[ d ]) + enVal_[ d ];
+ }
+ }
+ };
+
+
public:
virtual std::vector<double> computeTimeSteps () const
{
@@ -348,7 +459,7 @@ namespace Fem
return tmp;
}
- size_t leafElements() const
+ size_t numberOfElements() const
{
return elementCounter_;
}
@@ -485,17 +596,17 @@ namespace Fem
applyLocalImp( entity );
}
+ void setTroubledCellIndicator(TroubledCellIndicatorType indicator)
+ { // don't need a troubled cell indicator here
+ }
+
protected:
//! Perform the limitation on all elements.
void applyLocalImp(const EntityType& en) const
{
- //std::cout << "ScalingPass::applyLocalImp" << std::endl;
// timer for shock detection
Dune::Timer indiTime;
- // extract types
- enum { dim = EntityType :: dimension };
-
// check argument is not zero
assert( arg_ );
@@ -503,6 +614,9 @@ namespace Fem
// set entity to caller
caller().setEntity( en );
+ // bind function to entity
+ scaledFunction_.bind( en, spc_.order( en ) );
+
// get function to limit
const ArgumentFunctionType &U = *(std::get< argumentPosition >( *arg_ ));
@@ -512,7 +626,8 @@ namespace Fem
// get geometry
const Geometry& geo = en.geometry();
- RangeType enVal ;
+ // get reference to cell average
+ RangeType& enVal = scaledFunction_.enVal_;
bool limiter = false;
@@ -523,7 +638,7 @@ namespace Fem
limiter = true;
}
- RangeType theta( 1 );
+ RangeType& theta = scaledFunction_.resetTheta();
CornerPointSetType cornerquad( en );
if( ! checkPhysicalQuad( cornerquad, uEn, enVal, theta ) )
@@ -554,14 +669,23 @@ namespace Fem
limiter = true;
}
+ for (auto t : theta)
+ {
+ if (t<1) // if scaling would be done
+ limiter = true;
+ }
+
// scale function
if( limiter )
{
// get local funnction for limited values
DestLocalFunctionType limitEn = dest_->localFunction(en);
- // project scaled function to limitEn
- L2project(en, geo, enVal, uEn, theta, limitEn);
+ // set all dofs to zero
+ limitEn.clear();
+
+ const auto interpolation = spc_.interpolation( en );
+ interpolation( scaledFunction_, limitEn.localDofVector() );
// set indicator 1
discreteModel_.markIndicator();
@@ -583,20 +707,19 @@ namespace Fem
RangeType& theta ) const
{
// geometry and also use caching
- const int quadNop = quad.nop();
const EntityType& en = uEn.entity();
- tmpVal_.resize( quadNop );
+ auto& tmpVal = scaledFunction_.getValues( quad );
bool physical = true;
// evaluate uEn on all quadrature points
- uEn.evaluateQuadrature( quad , tmpVal_ );
+ uEn.evaluateQuadrature( quad , tmpVal );
- //const Geometry& geo = en.geometry();
+ const int quadNop = quad.nop();
for(int l=0; l<quadNop; ++l)
{
- const RangeType& u = tmpVal_[ l ];
+ const RangeType& u = tmpVal[ l ];
for( const auto& d : discreteModel_.model().limitedRange() )
{
@@ -609,8 +732,6 @@ namespace Fem
theta[ d ] = std::min( std::min( upper, lower ), theta[d ]);
}
- // warning!!! caching of geo can be more efficient
- //const DomainType& xgl = geo.global( quad[l] );
// check data
if ( ! discreteModel_.physical( en, quad.point(l), u ) )
{
@@ -623,78 +744,6 @@ namespace Fem
return physical;
}
- // check physicality on given quadrature
- template <class QuadratureType, class LocalFunctionImp>
- bool checkPhysicalQuad(const QuadratureType& quad,
- const LocalFunctionImp& uEn) const
- {
- // geometry and also use caching
- const int quadNop = quad.nop();
- const EntityType& en = uEn.entity();
-
- tmpVal_.resize( quadNop );
-
- // evaluate uEn on all quadrature points
- uEn.evaluateQuadrature( quad , tmpVal_ );
-
- //const Geometry& geo = en.geometry();
- for(int l=0; l<quadNop; ++l)
- {
- const RangeType& u = tmpVal_[ l ];
-
- // warning!!! caching of geo can be more efficient
- //const DomainType& xgl = geo.global( quad[l] );
- // check data
- if ( ! discreteModel_.physical( en, quad.point(l), u ) )
- {
- // notPhysical
- return false ;
- }
- }
- }
-
- //! check physicallity of data
- template <class LocalFunctionImp>
- bool checkPhysical(const EntityType& en,
- const Geometry& geo,
- const LocalFunctionImp& uEn) const
- {
- enum { dim = dimGrid };
- if( discreteModel_.hasPhysical() )
- {
-#if 1
- // use LagrangePointSet to evaluate on corners of the
- // geometry and also use caching
- return checkPhysicalQuad( CornerPointSetType( en ), uEn );
-#else
- {
- VolumeQuadratureType volQuad(en, volumeQuadOrd_ );
- if( ! checkPhysicalQuad(volQuad, uEn) ) return false;
- }
-
- const IntersectionIteratorType endnit = gridPart_.iend(en);
- for (IntersectionIteratorType nit = gridPart_.ibegin(en);
- nit != endnit; ++nit)
- {
- const IntersectionType& intersection = *nit;
- if( intersection.neighbor() && ! intersection.conforming() )
- {
- typedef typename FaceQuadratureType :: NonConformingQuadratureType NonConformingQuadratureType;
- NonConformingQuadratureType faceQuadInner(gridPart_,intersection, faceQuadOrd_, FaceQuadratureType::INSIDE);
- if( ! checkPhysicalQuad( faceQuadInner, uEn ) ) return false;
- }
- else
- {
- // conforming case
- FaceQuadratureType faceQuadInner(gridPart_,intersection, faceQuadOrd_, FaceQuadratureType::INSIDE);
- if( ! checkPhysicalQuad( faceQuadInner, uEn ) ) return false;
- }
- }
-#endif
- } // end physical
- return true;
- }
-
template <class LocalFunctionImp, class SpaceImp>
struct NumLinearBasis
{
@@ -732,82 +781,6 @@ namespace Fem
}
}
- // L2 projection
- template <class LocalFunctionImp>
- void L2project(const EntityType& en,
- const Geometry& geo,
- const RangeType& enVal,
- const LocalFunctionType& uEn,
- const RangeType& theta,
- LocalFunctionImp& limitEn ) const
- {
- enum { dim = dimGrid };
-
- // true if geometry mapping is affine
- const bool affineMapping = localMassMatrix_.affine();
-
- // set all dofs to zero
- limitEn.clear();
-
- // get quadrature for destination space order
- VolumeQuadratureType quad( en, volumeQuadratureOrder( en ) );
-
- //std::cout << globalMin << " " << globalMax << std::endl;
- //std::cout << minVal << " " << maxVal << std::endl;
-
- //old version that did not work well
- /*
- for( const auto& d : discreteModel_.model().limitedRange() )
- {
- double fst = 1.0;
- double sec = 1.0;
- if( std::abs( maxVal[ d ] - enVal[ d ] ) > 1e-10 )
- fst = (globalMax[ d ] - enVal[ d ])/(maxVal[ d ] - enVal[ d ]);
-
- if( std::abs( minVal[ d ] - enVal[ d ] ) > 1e-10 )
- sec = (globalMin[ d ] - enVal[ d ])/(minVal[ d ] - enVal[ d ]);
-
- theta[ d ] = std::min( std::min( std::abs( fst ), std::abs( sec ) ), double(1) );
- if( std::abs( theta[d ] ) < 1e-12 )
- theta[ d ] = 0;
- }
- */
-
- const int quadNop = tmpVal_.size();// quad.nop();
- assert( quadNop == int(quad.nop()) );
-
- // this has been done previously in checkPhysicalQuad
- // uEn.evaluateQuadrature( quad, tmpVal_ );
-
- for(int qP = 0; qP < quadNop ; ++qP)
- {
- RangeType &value = tmpVal_[ qP ];
-
- // \tilde{p}(x) = \theta (p(x) - \bar{u}) + \bar{u}
- // limitedRange contains all components that should be modified
- // default is 0,...,dimRange-1
- for( const auto& d : discreteModel_.model().limitedRange() )
- {
- value[ d ] = theta[ d ] * ( value[ d ] - enVal[ d ]) + enVal[ d ];
- }
-
- // get quadrature weight
- const double intel = (affineMapping) ?
- quad.weight(qP) : // affine case
- quad.weight(qP) * geo.integrationElement( quad.point(qP) ); // general case
-
- // quadrature scaling
- value *= intel;
- }
-
- // add all value to limitEn
- limitEn.axpyQuadrature( quad, tmpVal_ );
-
- // apply local inverse mass matrix for non-linear mappings
- if( !affineMapping )
- localMassMatrix_.applyInverse( en, limitEn );
- }
-
template <class BasisFunctionSetType, class PointType>
const RangeType& evaluateConstantBasis( const BasisFunctionSetType& basisSet,
const PointType& x ) const
@@ -839,7 +812,7 @@ namespace Fem
RangeType& val) const
{
bool notphysical = false;
- if( HierarchicalBasis< DiscreteFunctionSpaceType > :: v
+ if( Dune::Fem::Capabilities::isHierarchic< DiscreteFunctionSpaceType > :: v
&& localMassMatrix_.affine() )
{
// get point quadrature
@@ -928,6 +901,8 @@ namespace Fem
const DiscreteFunctionSpaceType& spc_;
GridPartType& gridPart_;
+ mutable ScaledFunction scaledFunction_;
+
const IndexSetType& indexSet_;
CornerPointSetContainerType cornerPointSetContainer_;
@@ -945,8 +920,6 @@ namespace Fem
mutable RangeType globalMin_;
mutable RangeType globalMax_ ;
- mutable std::vector< RangeType > tmpVal_ ;
-
mutable std::vector< RangeType > aver_ ;
// vector for stroing the information which elements have been computed already
diff --git a/dune/fem-dg/operator/limiter/smoothness.hh b/dune/fem-dg/operator/limiter/smoothness.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b4e44959f0f45bceb9220cdf84d43a8dd0965302
--- /dev/null
+++ b/dune/fem-dg/operator/limiter/smoothness.hh
@@ -0,0 +1,115 @@
+#include <cmath>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/dynvector.hh>
+#include <dune/common/dynmatrix.hh>
+#include <dune/fem/function/localfunction/const.hh>
+#include <dune/fem/space/shapefunctionset/orthonormal.hh>
+#include <dune/fem-dg/operator/limiter/indicatorbase.hh>
+
+template <class LF>
+double smoothnessIndicator(const LF& uLocal)
+{
+ using ONB = Dune::Fem::OrthonormalShapeFunctions<LF::dimDomain>;
+ const std::size_t R = LF::dimRange; // will be using density
+ const std::size_t P = uLocal.order();
+ assert(uLocal.size()/double(R) == ONB::size(P));
+ // actually some inverse mass matrix needed here
+ const double area = uLocal.entity().geometry().volume();
+ const double factor = 1/std::sqrt( area );
+ // 1D monoments vector
+ double q[P+1];
+ double b2[P+1];
+ double f = 0;
+ std::size_t k = ONB::size(0);
+ q[0] = uLocal[0]*uLocal[0]; // constant part not used
+ double l2norm2 = 0; // q[0];
+ for (std::size_t i=1;i<=P;++i)
+ {
+ q[i] = 0;
+ b2[i] = 0;
+ double nofMoments = ONB::size(i)-k;
+ assert(nofMoments == 1.);
+ for (;k<ONB::size(i);++k)
+ {
+ q[i] += uLocal[k*R]*uLocal[k*R] / nofMoments;
+ l2norm2 += uLocal[k*R]*uLocal[k*R] / nofMoments;
+ b2[i] += pow(1/double(i),2*P) / nofMoments;
+ f += pow(1/double(i),2*P) / nofMoments;
+ }
+ }
+#if 0
+ for (std::size_t i=0;i<=P;++i)
+ q[i] = std::sqrt( q[i] ) / factor;
+#endif
+ for (std::size_t i=1;i<=P;++i)
+ {
+ /*
+ std::cout << "====: " << i << " " << q[i] << " "
+ << l2norm2 << " " << b2[i]/f << " " << l2norm2*b2[i]/f
+ << " -> " << std::sqrt( q[i] + l2norm2*b2[i]/f ) / factor
+ << std::endl; */
+ q[i] = std::sqrt( q[i] + l2norm2*b2[i]/f ) / factor;
+ }
+ double maxQ = std::max( q[P], q[P-1] );
+ // find first 'significant' mode
+ std::size_t significant = 0;
+ for (std::size_t i=P;i>=1;--i)
+ {
+ maxQ = std::max(maxQ, q[i]);
+ if (maxQ>1e-14)
+ {
+ significant = i;
+ break;
+ }
+ }
+ if (significant==0) // constant, i.e., very smooth indeed
+ return 100;
+ if (significant==1) // not quite clear if this needs to be fixed -
+ // solution is linear and we do not have enough info to fit
+ return 10;
+
+ Dune::DynamicMatrix<double> matrix(significant,2);
+ Dune::DynamicVector<double> rhs(significant);
+ for (std::size_t r=significant; r-->0; )
+ {
+ maxQ = std::max(maxQ, q[r+1]);
+ rhs[r] = std::log( maxQ );
+ matrix[r][0] = 1;
+ matrix[r][1] = -std::log(double(r+1));
+ }
+ Dune::FieldMatrix<double,2,2> A;
+ Dune::FieldVector<double,2> b;
+ for (std::size_t r=0;r<2;++r)
+ {
+ for (std::size_t c=0;c<2;++c)
+ {
+ A[r][c] = 0;
+ for (std::size_t k=0;k<significant;++k)
+ A[r][c] += matrix[k][r]*matrix[k][c];
+ }
+ b[r] = 0;
+ for (std::size_t k=0;k<significant;++k)
+ b[r] += matrix[k][r]*rhs[k];
+ }
+ Dune::FieldVector<double,2> x;
+ A.solve(x,b);
+ double s = x[1];
+ assert(s==s);
+ // std::cout << " indicator= " << s << std::endl;
+ return s;
+}
+template <class DiscreteFunction>
+struct ModalSmoothnessIndicator
+: public Dune::Fem::TroubledCellIndicatorBase<DiscreteFunction>
+{
+ typedef typename DiscreteFunction :: LocalFunctionType LocalFunctionType;
+ double operator()( const DiscreteFunction& U, const LocalFunctionType& uEn) const override
+ {
+ double modalInd = smoothnessIndicator( uEn );
+ if( std::abs( modalInd ) > 1e-14 )
+ return 1.0 / modalInd ;
+ else
+ return 0.0;
+ }
+};
diff --git a/dune/fem-dg/pass/dginversemass.hh b/dune/fem-dg/pass/dginversemass.hh
index fcad09174fffe61041796a3a06f3cbebb4e1ade0..112838c274488e02d64010d262613ea466e50c7f 100644
--- a/dune/fem-dg/pass/dginversemass.hh
+++ b/dune/fem-dg/pass/dginversemass.hh
@@ -13,6 +13,8 @@
#include <dune/fem-dg/pass/pass.hh>
#include <dune/fem-dg/pass/discretemodel.hh>
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
+
namespace Dune
{
@@ -94,7 +96,7 @@ namespace Dune
explicit DGInverseMassPass ( PreviousPass &previousPass,
const DiscreteFunctionSpaceType &space )
: BaseType( previousPass, space, "DGInverseMassPass" ),
- localMassMatrix_( space, 2*space.order() ),
+ localMassMatrix_( space, [this](const int order) { return DefaultQuadrature<DiscreteFunctionSpaceType >::volumeOrder(order); } ),
arg_( 0 ),
dest_( 0 )
{}
@@ -106,7 +108,8 @@ namespace Dune
const int volQuadOrd = -1,
const int faceQuadOrd = -1)
: BaseType( previousPass, space, "DGInverseMassPass" ),
- localMassMatrix_( space, ( volQuadOrd == -1 ) ? (2*space.order()) : volQuadOrd ),
+ localMassMatrix_( space, [this,volQuadOrd](const int order) { return
+ ( volQuadOrd == -1 ) ? DefaultQuadrature<DiscreteFunctionSpaceType >::volumeOrder(order) : volQuadOrd; } ),
arg_( 0 ),
dest_( 0 )
{}
diff --git a/dune/fem-dg/pass/dgmasspass.hh b/dune/fem-dg/pass/dgmasspass.hh
index e405b602d69cfb7f6c36ce71bb47d43e3618798e..467867f7db73418406381876fc87dcec50d084fd 100644
--- a/dune/fem-dg/pass/dgmasspass.hh
+++ b/dune/fem-dg/pass/dgmasspass.hh
@@ -14,6 +14,7 @@
#include <dune/fem-dg/pass/dginversemass.hh>
#include <dune/fem-dg/misc/crs.hh>
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
namespace Dune
{
@@ -107,7 +108,8 @@ namespace Dune
: scalarSpace_( const_cast< GridPartType& > (key.gridPart()) ),
volumeQuadratureOrder_( 2 * scalarSpace_.order() ),
matrices_( Fem :: ThreadManager :: maxThreads() ),
- localMassMatrix_( scalarSpace_ ),
+ localMassMatrix_( scalarSpace_, [this](const int order) { return
+ DefaultQuadrature< ScalarDiscreteFunctionSpaceType >::volumeOrder(order); } ),
sequence_( -1 )
{
assert( Fem::ThreadManager::singleThreadMode() );
diff --git a/dune/fem-dg/pass/dgpass.hh b/dune/fem-dg/pass/dgpass.hh
index 0db741f100bedd1081496a4495a9375c97215c72..7a5e7fafc025ced5d39a19f230e619f159a7d4b7 100644
--- a/dune/fem-dg/pass/dgpass.hh
+++ b/dune/fem-dg/pass/dgpass.hh
@@ -19,6 +19,9 @@
#include <dune/fem-dg/pass/modelcaller.hh>
#include <dune/fem-dg/pass/discretemodel.hh>
#include <dune/fem/misc/compatibility.hh>
+#include <dune/fem-dg/operator/dg/defaultquadrature.hh>
+
+#include <dune/fem/io/parameter.hh>
namespace Dune
{
@@ -135,16 +138,34 @@ namespace Fem
PreviousPassType& pass,
const DiscreteFunctionSpaceType& spc,
const int volumeQuadOrd = -1,
- const int faceQuadOrd = -1 )
+ const int faceQuadOrd = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose() )
+ : LocalCDGPass( discreteModel, pass, spc, Dune::Fem::Parameter::container(), volumeQuadOrd, faceQuadOrd, verbose )
+ {}
+
+ //- Public methods
+ //! Constructor
+ //! \param discreteModel Actual discrete model definition (see dgdiscretemodels.hh)
+ //! \param pass Previous pass
+ //! \param spc Space belonging to the discrete function local to this pass
+ //! \param volumeQuadOrd defines the order of the volume quadrature which is by default 2* space polynomial order
+ //! \param faceQuadOrd defines the order of the face quadrature which is by default 2* space polynomial order
+ LocalCDGPass( DiscreteModelType& discreteModel,
+ PreviousPassType& pass,
+ const DiscreteFunctionSpaceType& spc,
+ const Dune::Fem::ParameterReader& /* parameter */,
+ const int volumeQuadOrd = -1,
+ const int faceQuadOrd = -1,
+ const bool verbose = Dune::Fem::Parameter::verbose() )
: BaseType(pass, spc),
caller_(),
discreteModel_(discreteModel),
- arg_(0),
- dest_(0),
+ arg_( nullptr ),
+ dest_( nullptr ),
spc_(spc),
gridPart_(spc_.gridPart()),
indexSet_(gridPart_.indexSet()),
- visited_(0),
+ visited_(),
updEn_(spc_),
updNeigh_(spc_),
valEnVec_( 20 ),
@@ -159,7 +180,7 @@ namespace Fem
#ifdef USE_CACHED_INVERSE_MASSMATRIX
localMassMatrix_( InverseMassProviderType :: getObject( MassKeyType( gridPart_ ) ) ),
#else
- localMassMatrix_( spc_ , volumeQuadOrd_ ),
+ localMassMatrix_( spc_ , [this](const int order) { return DefaultQuadrature<DiscreteFunctionSpaceType >::volumeOrder(order); } ),
#endif
reallyCompute_( true )
{
@@ -642,10 +663,10 @@ namespace Fem
// eval boundary Flux
wspeedS += caller().boundaryFlux(intersection,
- faceQuadInner,
- l,
- fluxEn,
- diffFluxEn)
+ faceQuadInner,
+ l,
+ fluxEn,
+ diffFluxEn)
* faceQuadInner.weight(l);
// apply weights
@@ -751,6 +772,7 @@ namespace Fem
flux *= intel;
}
+
// add values to local function
updEn.axpyQuadrature( volQuad, valJacEn_ );
}
@@ -921,13 +943,13 @@ namespace Fem
//! return default face quadrature order
static int defaultVolumeQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity ));
+ return DefaultQuadrature< DiscreteFunctionSpaceType >::volumeOrder( space.order( entity ) );
}
- //! return default face quadrature order
+ //! return default face quadrature
static int defaultFaceQuadratureOrder( const DiscreteFunctionSpaceType& space, const EntityType& entity )
{
- return (2 * space.order( entity )) + 1;
+ return DefaultQuadrature< DiscreteFunctionSpaceType >::faceOrder( space.order( entity ) );
}
protected:
diff --git a/dune/fem-dg/pass/threadpass.hh b/dune/fem-dg/pass/threadpass.hh
index 4f63f708cca26917b5e017949730fb288d8b3e0d..0aa7fa4c8492d43a7b2064aead4e8eec53ba18fa 100644
--- a/dune/fem-dg/pass/threadpass.hh
+++ b/dune/fem-dg/pass/threadpass.hh
@@ -213,6 +213,22 @@ namespace Fem
PreviousPassType& pass,
const DiscreteFunctionSpaceType& spc,
const int volumeQuadOrd = -1,
+ const int faceQuadOrd = -1)
+ : ThreadPass( discreteModel, pass, spc, Dune::Fem::Parameter::container(), volumeQuadOrd, faceQuadOrd )
+ {}
+
+ /** \brief Constructor
+ * \param discreteModel Actual discrete model
+ * \param pass Previous pass
+ * \param spc Space belonging to the discrete function local to this pass
+ * \param volumeQuadOrd defines the order of the volume quadrature which is by default 2* space polynomial order
+ * \param faceQuadOrd defines the order of the face quadrature which is by default 2* space polynomial order
+ */
+ ThreadPass(const DiscreteModelType& discreteModel,
+ PreviousPassType& pass,
+ const DiscreteFunctionSpaceType& spc,
+ const Dune::Fem::ParameterReader ¶meter,
+ const int volumeQuadOrd = -1,
const int faceQuadOrd = -1) :
BaseType(pass, spc),
delDofs_( spc ),
@@ -229,7 +245,8 @@ namespace Fem
faceQuadOrd_( faceQuadOrd ),
firstCall_( true ),
requireCommunication_( true ),
- sumComputeTime_( Fem :: Parameter :: getValue<bool>("fem.parallel.sumcomputetime", false ) )
+ parameter_( ¶meter ),
+ sumComputeTime_( parameter.getValue<bool>("fem.parallel.sumcomputetime", false ) )
{
// initialize quadratures before entering multithread mode
initializeQuadratures( spc, volumeQuadOrd, faceQuadOrd );
@@ -254,6 +271,17 @@ namespace Fem
#endif
// get information about communication
requireCommunication_ = passes_[ 0 ]->requireCommunication();
+ parameter_ = nullptr;
+ }
+
+ template <class TroubledCellIndicatorType>
+ void setTroubledCellIndicator( TroubledCellIndicatorType indicator )
+ {
+ const int maxThreads = Fem::ThreadManager::maxThreads();
+ for(int i=0; i<maxThreads; ++i)
+ {
+ passes_[ i ]->setTroubledCellIndicator( indicator );
+ }
}
static void initializeQuadratures( const DiscreteFunctionSpaceType& space, const int volQuadOrder, const int faceQuadOrder )
@@ -595,8 +623,17 @@ namespace Fem
{
// use separate discrete discrete model for each thread
discreteModels_[ thread ] = new DiscreteModelType( singleDiscreteModel_ );
- // create dg passes, the last bool disables communication in the pass itself
- passes_[ thread ] = new InnerPassType( *discreteModels_[ thread ], previousPass_, space(), volumeQuadOrd_, faceQuadOrd_ );
+ const bool verbose = Dune::Fem::Parameter::verbose() && Fem::ThreadManager::isMaster();
+ if( parameter_ )
+ {
+ // create dg passes, the last bool disables communication in the pass itself
+ passes_[ thread ] = new InnerPassType( *discreteModels_[ thread ], previousPass_, space(), *parameter_, volumeQuadOrd_, faceQuadOrd_, verbose );
+ }
+ else
+ {
+ // create dg passes, the last bool disables communication in the pass itself
+ passes_[ thread ] = new InnerPassType( *discreteModels_[ thread ], previousPass_, space(), volumeQuadOrd_, faceQuadOrd_, verbose );
+ }
return ;
}
@@ -732,6 +769,9 @@ namespace Fem
const int faceQuadOrd_;
mutable bool firstCall_;
bool requireCommunication_;
+
+ const Dune::Fem::ParameterReader* parameter_;
+
const bool sumComputeTime_;
};
//! @}
diff --git a/dune/fem-dg/python/operator.hh b/dune/fem-dg/python/operator.hh
index 9f0f46c15905bbc27a765d1100c0c38448cbd51a..cb66e247b3425155600b00c61d624c4b9a599cf1 100644
--- a/dune/fem-dg/python/operator.hh
+++ b/dune/fem-dg/python/operator.hh
@@ -73,7 +73,7 @@ namespace Dune
cls.def( "applyLimiter", []( Operator &self, DF &u) { self.applyLimiter(u); } );
// cls.def( "setTime", &Operator::setTime);
cls.def( "_setTime", &Operator::setTime);
- cls.def_property_readonly("timeStepEstimate",
+ cls.def_property_readonly("localTimeStepEstimate",
[](const Operator &self) -> std::tuple<double,double,double>
{ return {self.timeStepEstimate(),
self.explicitOperator().timeStepEstimate(),
diff --git a/dune/fem-dg/solver/rungekuttasolver.hh b/dune/fem-dg/solver/rungekuttasolver.hh
index 37fc487dad112b1c4373066028e6237374077118..b6261c0216f9f625e5ca31226221904a4dd296f0 100644
--- a/dune/fem-dg/solver/rungekuttasolver.hh
+++ b/dune/fem-dg/solver/rungekuttasolver.hh
@@ -291,7 +291,6 @@ namespace Fem
// create implicit or explicit ode solver
if( odeSolverType_ == 0 )
{
- std::cout << "Creating explicit ode solver" << std::endl;
solver = OdeSolversType :: createExplicitSolver( operator_, tp, rkSteps_, *param_, parameter, name_ );
}
else if (odeSolverType_ == 1)
@@ -346,7 +345,11 @@ namespace Fem
else
tpPtr_->init();
tpPtr_->init();
- std::cout << "cfl = " << double(tpPtr_->factor()) << ", T_0 = " << tpPtr_->time() << std::endl;
+
+ if( Dune::Fem::Parameter::verbose() )
+ {
+ std::cout << "cfl = " << double(tpPtr_->factor()) << ", T_0 = " << tpPtr_->time() << std::endl;
+ }
}
void initialize( const DestinationType& U )
diff --git a/dune/optim/common/smallobject.hh b/dune/optim/common/smallobject.hh
index f210d2fffb47c86f8370f07ed558d088c15d9c5f..c0b85477c1fb3743b565acdf6f88a5b3cc27e024 100644
--- a/dune/optim/common/smallobject.hh
+++ b/dune/optim/common/smallobject.hh
@@ -74,7 +74,7 @@ namespace Dune
static SmallObjectPool &instance ()
{
- static SmallObjectPool inst;
+ static thread_local SmallObjectPool inst;
return inst;
}
diff --git a/pydemo/euler/euler.py b/pydemo/euler/euler.py
index 93be7d81c6264bd923f389cb761c37fc6190227b..b70c3376ebc2a4bc98735822840803d4f06e2561 100644
--- a/pydemo/euler/euler.py
+++ b/pydemo/euler/euler.py
@@ -1,5 +1,6 @@
from ufl import *
from dune.common import FieldVector
+from dune.grid import reader
from dune.ufl import Space, GridFunction
from dune.fem.function import gridFunction
import numpy
@@ -128,7 +129,8 @@ def sod(dim=2,gamma=1.4):
#Model.domain=[[0, 0], [1, 0.25], [256, 64]]
#Model.domain=[[0, 0], [1, 0.25], [128, 32]]
Model.domain=[[0, 0], [1, 0.25], [64, 16]]
- Model.endTime=0.15
+ Model.domain=(reader.dgf, "grid2d_nonaffine.dgf")
+ Model.endTime=0.25
def chorin(gv,t):
gf = gv.function("chorin","chorin.hh", Vl,Vr,gamma,x0,t,name="chorin")
lgf = gf.localFunction() # this seems to fail?
@@ -137,7 +139,7 @@ def sod(dim=2,gamma=1.4):
lgf.bind(e)
return FieldVector( Model.toCons(lgf(x)) )
# return Model.toCons(gf(e,x))
- lf.plot()
+ # lf.plot()
return lf
Model.exact = lambda gv,t: chorin(gv,t)
Model.name="sod"
diff --git a/pydemo/euler/main.py b/pydemo/euler/main.py
index 6861621b6ca3c71208d6396cf615c37403f0fcdd..577bc3cd62d141409c577a2a26efd0c799315dee 100644
--- a/pydemo/euler/main.py
+++ b/pydemo/euler/main.py
@@ -1,6 +1,7 @@
-import mpi4py.rc
+#import mpi4py.rc
# mpi4py.rc.threaded = False
from dune.fem import parameter
+from dune.alugrid import aluCubeGrid
from dune.femdg.testing import run
# from euler import constant as problem
@@ -12,15 +13,13 @@ from euler import sod as problem
dim = 2
gamma = 1.4
-parameter.append({"fem.verboserank": -1})
+parameter.append({"fem.verboserank": 0})
primitive=lambda Model,uh: {"pressure": Model.toPrim(uh)[2]}
parameters = {"fem.ode.odesolver": "EX",
- "fem.timeprovider.factor": 0.25,
"fem.ode.order": 3,
- "femdg.limiter.admissiblefunctions": 1,
- "femdg.limiter.tolerance": 1,
- "femdg.limiter.epsilon": 1e-8}
+ "femdg.limiter.tolerance":1 }
+
#-----------------
# "dgadvectionflux.method": "EULER-HLLC", "EULER-HLL", "LLF"
# default value is 'LLF'
@@ -31,10 +30,16 @@ parameters = {"fem.ode.odesolver": "EX",
# 0 = only dg solution | 1 = only reconstruction | 2 = both
#-----------------
-run(problem(),
- startLevel=0, polOrder=2, limiter="default",
- primitive=primitive, saveStep=0.16, subsamp=2,
- dt=None,threading=False,grid="alucube", space="dgonb",
+Model = problem()
+Model.domain = aluCubeGrid(Model.domain,dimgrid=2)
+
+run(Model,
+ startLevel=3, polOrder=4, limiter="default",
+ primitive=primitive, saveStep=0.05, subsamp=2,
+ dt=None,threading=False,grid="alucube",
+ space="dglagrange",
+ # space=("dglagrange","lobatto"),
+ # space=("dglagrange","gauss"),
parameters=parameters)
# print(str(parameter))
diff --git a/pydemo/euler/testdg.py b/pydemo/euler/testdg.py
index ec0e7a639a3f5ecf8b22dd6717da87ccd0c7206a..d4c914d9a546cdd44db294aaf34b6f823f9f71cb 100644
--- a/pydemo/euler/testdg.py
+++ b/pydemo/euler/testdg.py
@@ -1,5 +1,7 @@
-import mpi4py.rc
+# import mpi4py.rc
# mpi4py.rc.threaded = False
+help("modules")
+
from dune.fem import parameter
from dune.femdg.testing import run
diff --git a/pydemo/euler/testfv0.py b/pydemo/euler/testfv0.py
index 04b78d05429677547f50f46a89da44b277162723..1befbcbb8d850f6542191f11049a0d2b202949d3 100644
--- a/pydemo/euler/testfv0.py
+++ b/pydemo/euler/testfv0.py
@@ -1,4 +1,4 @@
-import mpi4py.rc
+# import mpi4py.rc
# mpi4py.rc.threaded = False
from dune.fem import parameter
from dune.femdg.testing import run
diff --git a/pydemo/euler/testfv1.py b/pydemo/euler/testfv1.py
index 9b186a3043d2416037c3c6c9232cee32c5245f3c..b1123c08c32d1d246b242788ab012f40f6735b00 100644
--- a/pydemo/euler/testfv1.py
+++ b/pydemo/euler/testfv1.py
@@ -1,4 +1,4 @@
-import mpi4py.rc
+# import mpi4py.rc
# mpi4py.rc.threaded = False
from dune.fem import parameter
from dune.femdg.testing import run
diff --git a/pydemo/euler/testpolyfv1.py b/pydemo/euler/testpolyfv1.py
index b72cb5134f347554773018a2771492c28e24f8e6..ec88912ae4e536262b99502887decefc451fc506 100644
--- a/pydemo/euler/testpolyfv1.py
+++ b/pydemo/euler/testpolyfv1.py
@@ -1,4 +1,4 @@
-import mpi4py.rc
+# import mpi4py.rc
# mpi4py.rc.threaded = False
from dune.fem import parameter
from dune.femdg.testing import run
@@ -38,7 +38,7 @@ Model.exact = None
run(Model,
startLevel=0, limiter="default",
- primitive=primitive, saveStep=0.16, subsamp=2,
+ primitive=primitive, saveStep=0.16,
dt=None,threading=False,grid="polygon", space="finitevolume",
parameters=parameters)
diff --git a/python/dune/femdg/_operators.py b/python/dune/femdg/_operators.py
index fd06811049bcb742c7033789bf0fdfd8aa7c229a..b182b6c5706716d545a606a04b7179e468e97a8c 100644
--- a/python/dune/femdg/_operators.py
+++ b/python/dune/femdg/_operators.py
@@ -1,12 +1,15 @@
from __future__ import absolute_import, division, print_function, unicode_literals
-import sys
+import sys, hashlib
import logging
logger = logging.getLogger(__name__)
+from io import StringIO
from dune.common.checkconfiguration import assertHave, preprocessorAssert, ConfigurationError
+from dune.typeregistry import generateTypeName
from dune.generator import Constructor, Method
+from dune.generator.importclass import load as classLoad
from dune.common.hashit import hashIt
from dune.fem.operator import load
from dune.fem import parameter as parameterReader
@@ -83,31 +86,10 @@ def createOrderRedcution(domainSpace):
return load(includes, typeName, constructor).Operator( domainSpace )
#####################################################
-## fem-dg Operator
+## fem-dg models
#####################################################
from dune.femdg.patch import transform
-# create DG operator + solver (limiter = none,minmod,vanleer,superbee),
-# (diffusionScheme = cdg2,br2,ip,nipg,...)
-def femDGOperator(Model, space,
- limiter="minmod",
- advectionFlux="default",
- diffusionScheme = "cdg2", threading=False,
- initialTime=0.0, parameters=None):
- virtualize = False
- import dune.create as create
-
- includes = []
-
- if limiter is None or limiter is False:
- limiter = "unlimited"
-
- if limiter.lower() == "default":
- limiter = "minmod"
-
- # TODO: does this make sense - if there is no diffusion then it doesn't
- # matter and with diffusion using 'none' seems a bad idea?
- if diffusionScheme is None or diffusionScheme is False:
- diffusionScheme = "none"
+def femDGModels(Model, space, initialTime=0):
u = TrialFunction(space)
v = TestFunction(space)
@@ -144,13 +126,124 @@ def femDGOperator(Model, space,
diffModel += inner(as_vector(Model.S_impl(t,x,u,grad(u))),v)*dx
print("Model.S_s is deprecated. Use S_i instead!")
- advModel = create.model("elliptic",space.grid, advModel,
+ from dune.fem.model import elliptic
+ virtualize = False
+
+ advModel = elliptic(space.grid, advModel,
modelPatch=transform(Model,space,t,"Adv"),
virtualize=virtualize)
- diffModel = create.model("elliptic",space.grid, diffModel,
+ diffModel = elliptic(space.grid, diffModel,
modelPatch=transform(Model,space,t,"Diff"),
virtualize=virtualize)
+ Model._ufl = {"u":u,"v":v,"n":n,"x":x,"t":t}
+
+ return [Model,advModel,diffModel]
+
+#####################################################
+## fem-dg Operator
+#####################################################
+# create DG operator + solver (limiter = none,default,minmod,vanleer,superbee,lp,scaling),
+# (diffusionScheme = cdg2,cdg,br2,ip,nipg,bo)
+def femDGOperator(Model, space,
+ limiter="default",
+ advectionFlux="default",
+ diffusionScheme = "cdg2",
+ threading=False,
+ defaultQuadrature=True,
+ codegen=True,
+ initialTime=0.0, parameters=None):
+ """ create DG operator + ODE solver
+
+ Args:
+ Model: analytical model describing the PDE and auxiliary functionality
+ space: discrete function space (DG space)
+ limiter: choice of limiter stabilization, possible values are
+ none,default,minmod,vanleer,superbee,lp,scaling
+ advectionFlux: choice of numerical flux for advective parts
+ default is local Lax-Friedrichs
+ diffusionScheme: choice of numerical flux for diffusive parts
+ possible choices are cdg2(default),cdg,br2,ip,nipg,bo
+ threading: enable shared memory parallelization
+ defaultQuadrature: use quadratures that generically fit to the space
+ codegen: enable optimized code for evaluation and interpolation
+ initialTime: T_0, default is 0.0
+ parameters: Additional parameter passed to the DG operator, limiter and ODE solvers
+ Returns:
+ DGOperator
+ """
+
+ includes = []
+
+ # obtain information about limiter interface before Model
+ # is augmented with default implementations
+ hasScalingInterface = hasattr(Model,"lowerBound") or hasattr(Model,"upperBound") or hasattr(Model,"physical")
+ hasLimiterInterface = (hasattr(Model,"jump") and hasattr(Model,"velocity")) or hasattr(Model,"physical")
+
+ if type(Model)==list or type(Model)==tuple:
+ advModel = Model[1]
+ diffModel = Model[2]
+ Model = Model[0]
+ else:
+ Model, advModel, diffModel = femDGModels(Model,space,initialTime)
+
+ hasAdvFlux = hasattr(Model,"F_c")
+ hasDiffFlux = hasattr(Model,"F_v")
+ hasStiffSource = hasattr(Model,"S_i") or hasattr(Model,"S_s")
+ hasNonStiffSource = hasattr(Model,"S_e") or hasattr(Model,"S_ns")
+
+ virtualize = False
+
+ defaultLimiter = limiter == "default"
+
+ if limiter is None or limiter is False:
+ limiter = "unlimited"
+
+ if type(limiter) in [list,tuple]:
+ limiterIndicator = limiter[1]
+ limiter = limiter[0]
+ else:
+ limiterIndicator = None
+
+ limiterstr = limiter
+ if limiter.lower() == "default":
+ # check for limiter interface implementation
+ if not hasLimiterInterface:
+ print("\nfemDGOperator: Limiter selected but limiter interface (jump,velocity,physical) missing in Model!", flush=True)
+ print("femDGOperator: Falling back to unlimited!\n", flush=True)
+ limiter = "unlimited"
+ else:
+ # default is minmod which can be either lp-minmod or muscl-minmod
+ limiter = "minmod"
+ limiterstr = limiter if space.grid.type.isSimplex else "lp"
+ # force default values for how reconstruction is done
+ parameters["femdg.limiter.admissiblefunctions"] = "default"
+
+ if limiter.lower() == "lp":
+ limiter = "minmod"
+ # force default values for how reconstruction is done
+ parameters["femdg.limiter.admissiblefunctions"] = "lp"
+
+ if limiter.lower() == "scaling":
+ # check for scaling limiter interface
+ if not hasScalingInterface:
+ raise KeyError(\
+ "femDGOperator: ScalingLimiter selected but scaling limiter interface (lowerBound,upperBound,physical) missing in Model!\n")
+ elif limiter.lower() != "unlimited":
+ # check for scaling limiter interface
+ if not hasScalingInterface:
+ raise KeyError(\
+ "femDGOperator: MUSCL type stabilization selected but limiter interface (jump,velocity,physical) missing in Model!\n")
+
+ if space.grid.comm.rank == 0:
+ limiterstr = "default(" + limiterstr + ")" if defaultLimiter else limiterstr
+ print("femDGOperator: Limiter =",limiterstr)
+
+ # TODO: does this make sense - if there is no diffusion then it doesn't
+ # matter and with diffusion using 'none' seems a bad idea?
+ if diffusionScheme is None or diffusionScheme is False:
+ diffusionScheme = "none"
+
spaceType = space._typeName
if virtualize:
@@ -197,13 +290,13 @@ def femDGOperator(Model, space,
# at the moment this is always the same type (depending on
# model._typeName) so could be done by only providing the header
# file in the dg operator construction method
- from dune.typeregistry import generateTypeName
clsName,includes = generateTypeName("Dune::Fem::DGAdvectionFlux",advModel._typeName,"Dune::Fem::AdvectionFlux::Enum::userdefined")
advectionFlux = advectionFlux(advModel,clsName,includes)
includes += advectionFlux._includes
- #elif advectionFlux.lower().find(".h") >= 0:
- # advFluxId = "Dune::Fem::AdvectionFlux::Enum::userdefined"
- # includes += [ advectionFlux ]
+ elif hasattr(advectionFlux,"_typeName"):
+ advFluxId = "Dune::Fem::AdvectionFlux::Enum::userdefined"
+ advectionFluxIsCallable = True
+ includes += advectionFlux._includes
else:
# if dgadvectionflux.method has been selected, then use general flux,
# otherwise default to LLF flux
@@ -242,7 +335,8 @@ def femDGOperator(Model, space,
elif limiter.lower() != "minmod":
raise ValueError("limiter "+limiter+" not recognised")
- signature = (advFluxId,diffusionScheme,threading,solverId,formId,limiterId,limiterFctId,advFluxId,diffFluxId,)
+ signature = (advFluxId,diffusionScheme,threading,solverId,formId,
+ limiterId,limiterFctId,advFluxId,diffFluxId,str(defaultQuadrature))
additionalClass = "Additional_"+hashIt(str(signature))
struct = Struct(additionalClass, targs=['class FunctionSpace'])
struct.append(TypeAlias('DomainType','typename FunctionSpace::DomainType'))
@@ -258,8 +352,7 @@ def femDGOperator(Model, space,
limitedDimRange = "FunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len( limiterModifiedDict.items() )
limitedDimRange = str(count)
struct.append([Declaration(
Variable("const int", "limitedDimRange = " + limitedDimRange),
@@ -303,6 +396,9 @@ def femDGOperator(Model, space,
struct.append([Declaration(
Variable("const Dune::Fem::DiffusionFlux::Enum", "diffFluxId = " + diffFluxId),
static=True)])
+ struct.append([Declaration(
+ Variable("const bool", "defaultQuadrature"), initializer=defaultQuadrature,
+ static=True)])
writer = SourceWriter(StringWriter())
writer.emit([struct])
@@ -329,33 +425,50 @@ def femDGOperator(Model, space,
_, domainFunctionIncludes, domainFunctionType, _, _, _ = space.storage
base = 'Dune::Fem::SpaceOperatorInterface< ' + domainFunctionType + '>'
- estimateMark = Method('estimateMark', '''[]( DuneType &self, const typename DuneType::DestinationType &u, const double dt) { self.estimateMark(u, dt); }''' );
+ estimateMark = Method('estimateMark', '''[]( DuneType &self, const typename DuneType::DestinationType &u, const double dt) { self.estimateMark(u, dt); }''' )
+
+ includes += ["dune/fem-dg/operator/limiter/indicatorbase.hh"]
+ setIndicator = Method('setTroubledCellIndicator',
+ args=['DuneType &self, typename DuneType::TroubledCellIndicatorType indicator'],
+ body=['self.setTroubledCellIndicator(indicator);'],
+ extra=['pybind11::keep_alive<0,1>()'])
+ gridSizeInterior = Method('gridSizeInterior', '&DuneType::gridSizeInterior')
+
+ order = space.order
+ if codegen:
+ codegen = [space,range(2,2*order+2),range(2,2*order+2)]
+ else:
+ codegen = None
if parameters is not None:
if advectionFluxIsCallable:
- op = load(includes, typeName, estimateMark,
- baseClasses = [base],
+ op = load(includes, typeName, estimateMark, setIndicator, gridSizeInterior,
+ baseClasses=[base],
+ codegen=codegen,
preamble=writer.writer.getvalue()).\
Operator( space, advModel, diffModel, advectionFlux, parameters=parameters )
else:
- op = load(includes, typeName, estimateMark,
+ op = load(includes, typeName, estimateMark, setIndicator, gridSizeInterior,
baseClasses = [base],
+ codegen=codegen,
preamble=writer.writer.getvalue()).\
Operator( space, advModel, diffModel, parameters=parameters )
else:
if advectionFluxIsCallable:
- op = load(includes, typeName, estimateMark,
+ op = load(includes, typeName, estimateMark, setIndicator, gridSizeInterior,
baseClasses = [base],
+ codegen=codegen,
preamble=writer.writer.getvalue()).\
Operator( space, advModel, diffModel, advectionFlux )
else:
- op = load(includes, typeName, estimateMark,
+ op = load(includes, typeName, estimateMark, setIndicator, gridSizeInterior,
baseClasses = [base],
+ codegen=codegen,
preamble=writer.writer.getvalue()).\
Operator( space, advModel, diffModel )
- op._t = t
- op.time = t.value
+ op._t = Model._ufl["t"]
+ op.time = Model._ufl["t"].value
op.models = [advModel,diffModel]
op.space = space
def setTime(self,time):
@@ -371,8 +484,40 @@ def femDGOperator(Model, space,
op.stepTime = stepTime.__get__(op)
op._hasAdvFlux = hasAdvFlux
op._hasDiffFlux = hasDiffFlux
+ if limiterIndicator is not None:
+ op.setTroubledCellIndicator(limiterIndicator)
return op
+def smoothnessIndicator(clsName, includes, u_h, ctorArgs=None):
+ if ctorArgs is None: ctorArgs=[]
+ baseName,_ = generateTypeName("Dune::Fem::TroubledCellIndicatorBase",u_h)
+ return classLoad(clsName, includes,*ctorArgs, baseClasses=[baseName],
+ holder="std::shared_ptr")
+def advectionNumericalFlux(clsName ,includes, advModel, additionalArgs=None):
+ if additionalArgs is None: additionalArgs=[]
+ code = '''
+namespace Dune
+{
+ namespace Fem
+ {
+ template <>
+ struct DGAdvectionFlux< XXXADV_MODELXXX, AdvectionFlux::Enum::userdefined >
+ : public XXXIMPL_MODELXXX
+ {
+ template< class ... Args>
+ DGAdvectionFlux( Args&&... args )
+ : XXXIMPL_MODELXXX( std::forward<Args>(args)... ) {}
+ std::string name() const {return "user defined flux";}
+ };
+ }
+}
+'''
+ code = code.replace("XXXADV_MODELXXX",advModel._typeName)
+ code = code.replace("XXXIMPL_MODELXXX",clsName)
+ clsName,includesA = generateTypeName("Dune::Fem::DGAdvectionFlux",advModel,
+ "Dune::Fem::AdvectionFlux::Enum::userdefined")
+ return classLoad(clsName, includes+includesA+[StringIO(code)], advModel, *additionalArgs)
+
# RungeKutta solvers
def rungeKuttaSolver( fullOperator, imex='EX', butchertable=None, parameters={} ):
@@ -410,9 +555,11 @@ def rungeKuttaSolver( fullOperator, imex='EX', butchertable=None, parameters={}
setTimeStepSize = Method('setTimeStepSize', '&DuneType::setTimeStepSize')
deltaT = Method('deltaT', '&DuneType::deltaT')
- return load(includes, typeName, constructor, solve, setTimeStepSize, deltaT).Operator(
- fullOperator.fullOperator,
- fullOperator.explicitOperator,
- fullOperator.implicitOperator,
- imexId,
- parameters=parameters)
+ return load(includes, typeName,
+ constructor, solve, setTimeStepSize, deltaT,
+ codegen=fullOperator.codegen).\
+ Operator( fullOperator.fullOperator,
+ fullOperator.explicitOperator,
+ fullOperator.implicitOperator,
+ imexId,
+ parameters=parameters )
diff --git a/python/dune/femdg/patch.py b/python/dune/femdg/patch.py
index 599091f293d8a23cb2687ed494f0a05a2ffea78f..65d132fb0efd5ee1feb963da3d05138ac0e1a38d 100644
--- a/python/dune/femdg/patch.py
+++ b/python/dune/femdg/patch.py
@@ -1,9 +1,12 @@
-from __future__ import division, print_function, unicode_literals
-# from dune.ufl.codegen import generateMethod
-from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal, Coefficient, replace, diff, as_vector
+from functools import reduce
+
+from ufl import grad, TrialFunction, SpatialCoordinate, FacetNormal,\
+ Coefficient, replace, diff, as_vector,\
+ conditional
from ufl.core.expr import Expr
-from dune.source.cplusplus import Variable, UnformattedExpression, AccessModifier, Declaration
+from dune.source.cplusplus import Variable, UnformattedExpression,\
+ AccessModifier, Declaration, Method
from ufl.algorithms import expand_compounds, expand_derivatives, expand_indices, expand_derivatives
def uflExpr(Model,space,t):
@@ -11,6 +14,30 @@ def uflExpr(Model,space,t):
n = FacetNormal(space.cell())
x = SpatialCoordinate(space.cell())
+ physicalBound = None
+ lowerBound = getattr(Model,"lowerBound",None)
+ if lowerBound is None:
+ lowerBound = space.dimRange*["std::numeric_limits<double>::min()"]
+ else:
+ lowerCond = [conditional(u[i]>=lowerBound[i],1,0)
+ for i in range(len(lowerBound))
+ if lowerBound[i] is not None ]
+ if lowerCond != []:
+ physicalBound = reduce( (lambda x,y: x*y), lowerCond )
+ else: lowerBound=space.dimRange*["std::numeric_limits<double>::min()"]
+ upperBound = getattr(Model,"upperBound",None)
+ if upperBound is None:
+ upperBound = space.dimRange*["std::numeric_limits<double>::max()"]
+ else:
+ upperCond = [conditional(u[i]<=upperBound[i],1,0)
+ for i in range(len(upperBound))
+ if upperBound[i] is not None ]
+ if upperCond != []:
+ physicalBound_ = reduce( (lambda x,y: x*y), upperCond )
+ physicalBound = physicalBound_ if physicalBound is None else\
+ physicalBound*physicalBound_
+ else: upperBound=space.dimRange*["std::numeric_limits<double>::min()"]
+
maxSpeed = getattr(Model,"maxLambda",None)
if maxSpeed is not None:
maxSpeed = maxSpeed(t,x,u,n)
@@ -20,9 +47,14 @@ def uflExpr(Model,space,t):
diffusionTimeStep = getattr(Model,"maxDiffusion",None)
if diffusionTimeStep is not None:
diffusionTimeStep = diffusionTimeStep(t,x,u)
- physical = getattr(Model,"physical",None)
- if physical is not None:
+ physical = getattr(Model,"physical",True)
+ if not isinstance(physical,bool):
physical = physical(t,x,u)
+ if physicalBound is not None:
+ physical = physical*physicalBound
+ else:
+ if physicalBound is not None:
+ physical = physical*physicalBound
# TODO: jump is problematic since the coefficient 'w' causes issues -
# no predefined when extracting required Coefficients so needs fixing
# So jump is not returned by this method and is constructed again in
@@ -61,7 +93,7 @@ def uflExpr(Model,space,t):
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"
+ 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)
@@ -78,12 +110,13 @@ def uflExpr(Model,space,t):
limitedDimRange = "DFunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len(limiterModifiedDict.items())
+ # for id,f in limiterModifiedDict.items(): count += 1
limitedDimRange = str(count)
# jump = None # TODO: see comment above
return maxSpeed, velocity, diffusionTimeStep, physical, jump,\
- boundaryAFlux, boundaryDFlux, boundaryValue, hasBoundaryValue
+ boundaryAFlux, boundaryDFlux, boundaryValue, hasBoundaryValue,\
+ physicalBound
def codeFemDg(self):
code = self._code()
@@ -100,6 +133,30 @@ def codeFemDg(self):
n = FacetNormal(space.cell())
x = SpatialCoordinate(space.cell())
+ physicalBound = None
+ lowerBound = getattr(Model,"lowerBound",None)
+ if lowerBound is None:
+ lowerBound = space.dimRange*["std::numeric_limits<double>::min()"]
+ else:
+ lowerCond = [conditional(u[i]>=lowerBound[i],1,0)
+ for i in range(len(lowerBound))
+ if lowerBound[i] is not None ]
+ if lowerCond != []:
+ physicalBound = reduce( (lambda x,y: x*y), lowerCond )
+ else: lowerBound=space.dimRange*["std::numeric_limits<double>::min()"]
+ upperBound = getattr(Model,"upperBound",None)
+ if upperBound is None:
+ upperBound = space.dimRange*["std::numeric_limits<double>::max()"]
+ else:
+ upperCond = [conditional(u[i]<=upperBound[i],1,0)
+ for i in range(len(upperBound))
+ if upperBound[i] is not None ]
+ if upperCond != []:
+ physicalBound_ = reduce( (lambda x,y: x*y), upperCond )
+ physicalBound = physicalBound_ if physicalBound is None else\
+ physicalBound*physicalBound_
+ else: upperBound=space.dimRange*["std::numeric_limits<double>::min()"]
+
# TODO come up with something better!
hasGamma = getattr(Model,"gamma",None)
code.append([Declaration(
@@ -158,10 +215,15 @@ def codeFemDg(self):
targs=['class Entity, class Point, class T'], const=True,
predefined=predefined)
- # QUESTION: should `physical` actually depend on x? Perhaps even t?
physical = getattr(Model,"physical",True)
if not isinstance(physical,bool):
physical = physical(t,x,u)
+ if physicalBound is not None:
+ physical = physical*physicalBound
+ else:
+ if physicalBound is not None:
+ physical = physical*physicalBound
+
self.generateMethod(code, physical,
'double', 'physical',
args=['const Entity &entity', 'const Point &x',
@@ -285,15 +347,24 @@ def codeFemDg(self):
limitedDimRange = "DFunctionSpace :: dimRange"
else:
limiterModified = {}
- count = 0
- for id,f in limiterModifiedDict.items(): count += 1
+ count = len(limiterModifiedDict.items())
+ # for id,f in limiterModifiedDict.items(): count += 1
limitedDimRange = str(count)
self.generateMethod(code, limiterModified,
'void', 'limitedRange',
args=['LimitedRange& limRange'],
targs=['class LimitedRange'], const=True, evalSwitch=False,
predefined=predefined)
-
+ obtainBounds = Method("void","obtainBounds",
+ targs=['class LimitedRange'],
+ args=["LimitedRange& globalMin","LimitedRange& globalMax"],const=True)
+ obtainBounds.append("globalMin = {"
+ + ', '.join(map(str, lowerBound))
+ +"};")
+ obtainBounds.append("globalMax = {"
+ + ', '.join(map(str, upperBound))
+ +"};")
+ code.append( [obtainBounds] )
return code
def transform(Model,space,t, name="" ):
diff --git a/python/dune/femdg/rk.py b/python/dune/femdg/rk.py
index 4db0e984509e393761d192e0277f79fae8ca6a32..9514364b10d2a1678f27a181d4cf106f49145e28 100644
--- a/python/dune/femdg/rk.py
+++ b/python/dune/femdg/rk.py
@@ -12,11 +12,17 @@ class FemDGStepper:
self.rkScheme.setTimeStepSize(dt)
self.rkScheme.solve(u)
return self.rkScheme.deltaT()
-def femdgStepper(*,order=None,rkType=None,operator=None,cfl=None,parameters=True):
- if parameters is True:
- parameters = {}
- elif parameters is False:
- parameters = None
+ @property
+ def deltaT(self):
+ return self.rkScheme.deltaT()
+ @deltaT.setter
+ def deltaT(self,dt):
+ self.rkScheme.setTimeStepSize(dt)
+
+def femdgStepper(*,order=None,rkType=None,operator=None,cfl=0.45,parameters=True):
+ if parameters is True: parameters = {}
+ elif parameters is False: parameters = None
+ if rkType == "default": rkType = None
def _femdgStepper(op,cfl=None):
if parameters is not None:
if not "fem.timeprovider.factor" in parameters:
@@ -148,7 +154,7 @@ class RungeKutta:
self.op.stepTime(0,0)
self.op(u,self.k[0])
if dt is None and self.dt is None:
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -158,12 +164,12 @@ class RungeKutta:
self.tmp.axpy(dt*self.A[i][j],self.k[j])
self.op.stepTime(self.c[i],dt)
self.op(self.tmp,self.k[i])
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
else:
if dt is None and self.dt is None:
self.op.stepTime(0,0)
self.op(u,self.k[0])
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -173,7 +179,7 @@ class RungeKutta:
self.tmp.axpy(dt*self.A[i][j],self.k[j])
self.op.stepTime(self.c[i],dt)
self.op(self.tmp,self.k[i]) # this seems like a good initial guess for dt small
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
self.helmholtz.alpha = dt*self.A[i][i]
self.helmholtz.solve(baru=self.tmp,target=self.k[i])
@@ -181,7 +187,7 @@ class RungeKutta:
u.axpy(dt*self.b[i],self.k[i])
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
class Heun(RungeKutta):
def __init__(self, op, cfl=None):
A = [[0,0],
@@ -225,7 +231,7 @@ class ImplSSP2: # with stages=1 same as above - increasing stages does not impro
if dt is None and self.dt is None:
self.op.stepTime(0,0)
self.op(u, self.tmp)
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -237,18 +243,18 @@ class ImplSSP2: # with stages=1 same as above - increasing stages does not impro
for i in range(2,self.stages+1):
self.op.stepTime(self.c(i),dt)
self.op(self.q2, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
self.q2.axpy(dt*self.mu21, self.tmp)
self.q2.assgin(tmp)
self.helmholtz.solve(self.tmp, self.q2)
u.as_numpy[:] *= (1-self.lamsps)
u.axpy(self.lamsps, self.q2)
self.op(self.q2, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(dt*self.musps, self.tmp)
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
class ExplSSP2:
def __init__(self,stages,op,cfl=None):
self.op = op
@@ -263,7 +269,7 @@ class ExplSSP2:
if dt is None and self.dt is None:
self.op.stepTime(0,0)
self.op(u, self.tmp)
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -272,17 +278,17 @@ class ExplSSP2:
for i in range(1,self.stages):
self.op.stepTime(self.c(i),dt)
self.op(u,self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(fac, self.tmp)
self.op.stepTime(self.c(i),dt)
self.op(u,self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.as_numpy[:] *= (self.stages-1)/self.stages
u.axpy(dt/self.stages, self.tmp)
u.axpy(1/self.stages, self.q2)
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
def ssp2(stages,explicit=True):
if explicit:
return lambda op,cfl=None: ExplSSP2(stages,op,cfl)
@@ -313,7 +319,7 @@ class ExplSSP3:
if dt is None and self.dt is None:
self.op.stepTime(0,0)
self.op(u, self.tmp)
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -322,14 +328,14 @@ class ExplSSP3:
while i <= (self.n-1)*(self.n-2)/2:
self.op.stepTime(self.c(i),dt)
self.op(u,self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(fac, self.tmp)
i += 1
self.q2.assign(u)
while i <= self.n*(self.n+1)/2:
self.op.stepTime(self.c(i),dt)
self.op(u,self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(fac, self.tmp)
i += 1
u.as_numpy[:] *= (self.n-1)/(2*self.n-1)
@@ -337,12 +343,12 @@ class ExplSSP3:
while i <= self.stages:
self.op.stepTime(self.c(i),dt)
self.op(u,self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(fac, self.tmp)
i += 1
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
class ImplSSP3:
def __init__(self,stages,op,cfl=None):
self.stages = stages
@@ -372,7 +378,7 @@ class ImplSSP3:
if dt is None and self.dt is None:
# self.op.stepTime(0,0)
self.op(u, self.tmp)
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -383,18 +389,18 @@ class ImplSSP3:
for i in range(2,self.stages+1):
# self.op.stepTime(self.c(i),dt)
self.op(self.q2, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
self.q2.axpy(dt*self.mu21, self.tmp)
self.helmholtz.solve(self.q2,self.tmp)
self.q2.assign(self.tmp)
u.as_numpy[:] *= (1-self.lamsps)
u.axpy(self.lamsps, self.q2)
self.op(self.q2, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(dt*self.musps, self.tmp)
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
def ssp3(stages,explicit=True):
if explicit:
return lambda op,cfl=None: ExplSSP3(stages,op,cfl)
@@ -415,7 +421,7 @@ class ExplSSP4_10:
if dt is None and self.dt is None:
self.op.stepTime(0,0)
self.op(u, self.tmp)
- dt = self.op.timeStepEstimate[0]*self.cfl
+ dt = self.op.localTimeStepEstimate[0]*self.cfl
elif dt is None:
dt = self.dt
self.dt = 1e10
@@ -425,7 +431,7 @@ class ExplSSP4_10:
while i <= 5:
self.op.stepTime(self.c(i), dt)
self.op(u, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(dt/6, self.tmp)
i += 1
@@ -437,16 +443,16 @@ class ExplSSP4_10:
while i <= 9:
self.op.stepTime(self.c(i), dt)
self.op(u, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.axpy(dt/6, self.tmp)
i += 1
self.op.stepTime(self.c(i), dt)
self.op(u, self.tmp)
- self.dt = min(self.dt, self.op.timeStepEstimate[0]*self.cfl)
+ self.dt = min(self.dt, self.op.localTimeStepEstimate[0]*self.cfl)
u.as_numpy[:] *= 3/5
u.axpy(1, self.q2)
u.axpy(dt/10, self.tmp)
self.op.applyLimiter( u )
self.op.stepTime(0,0)
- return dt
+ return self.op.space.grid.comm.min(dt)
diff --git a/python/dune/femdg/testing.py b/python/dune/femdg/testing.py
index 69f63f4b42a0d18dcba8f13530456e5bc3f59b0e..1c85ebb191c0a2126cdc22510190ab39904b669b 100644
--- a/python/dune/femdg/testing.py
+++ b/python/dune/femdg/testing.py
@@ -21,6 +21,11 @@ def run(Model, Stepper=None,
exact = Model.exact
except:
exact = None
+
+ if grid == "polygon":
+ subsamp=None
+ print('Disabling subsampling for PolygonGrid')
+
print("*************************************")
print("**** Running simulation",Model.name)
print("*************************************")
@@ -58,7 +63,10 @@ def run(Model, Stepper=None,
# create discrete function space
try:
- if space.lower() == "finitevolume":
+ if type(space) in [list,tuple]:
+ space = create.space( space[0], grid, dimRange=dimR, order=polOrder,
+ pointType=space[1] )
+ elif space.lower() == "finitevolume":
space = create.space( space, grid, dimRange=dimR)
else:
space = create.space( space, grid, order=polOrder, dimRange=dimR)
@@ -67,7 +75,10 @@ def run(Model, Stepper=None,
if modifyModel is not None:
Model = modifyModel(Model,space)
- operator = femDGOperator(Model, space, limiter=limiter, threading=threading, parameters=parameters )
+ operator = femDGOperator(Model, space,
+ limiter=limiter, threading=threading, parameters=parameters,
+ defaultQuadrature=True
+ )
stepper = Stepper(operator, cfl)
# create and initialize solution
u_h = space.interpolate(Model.initial, name='u_h')
@@ -120,7 +131,7 @@ def run(Model, Stepper=None,
sys.exit(1)
if t > saveTime:
print('[',tcount,']','dt = ', dt, 'time = ',t,
- 'dtEst = ',operator.timeStepEstimate,
+ 'dtEst = ',operator.localTimeStepEstimate,
'elements = ',grid.size(0), flush=True )
try:
velo0.setConstant("time",[t])