From b1b142dd0450cb376f446839978196cf6c1772be Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ren=C3=A9=20He=C3=9F?= <rene.hess@iwr.uni-heidelberg.de>
Date: Thu, 16 Nov 2017 15:32:06 +0100
Subject: [PATCH] Fix driver generation for taylor green vortex
---
python/dune/perftool/compile.py | 6 ++-
python/dune/perftool/options.py | 5 ++-
.../dune/perftool/pdelab/driver/__init__.py | 14 +++++++
.../perftool/pdelab/driver/constraints.py | 12 ++++--
python/dune/perftool/pdelab/driver/error.py | 24 +++++++++---
.../pdelab/driver/gridfunctionspace.py | 25 ++++++++++--
.../perftool/pdelab/driver/instationary.py | 28 ++++++++-----
python/dune/perftool/pdelab/driver/solve.py | 21 ++++++++--
test/navier-stokes/howto/taylor-green.cc | 33 ++++++++--------
test/navier-stokes/howto/taylor-green.ini | 2 +-
.../navierstokes_2d_dg_quadrilateral.mini | 13 +++++--
.../navierstokes_2d_dg_quadrilateral.ufl | 39 +++++++++----------
12 files changed, 156 insertions(+), 66 deletions(-)
diff --git a/python/dune/perftool/compile.py b/python/dune/perftool/compile.py
index 101ed2a7..353eda86 100644
--- a/python/dune/perftool/compile.py
+++ b/python/dune/perftool/compile.py
@@ -16,7 +16,10 @@ from ufl.algorithms.formfiles import interpret_ufl_namespace
from dune.perftool.generation import (delete_cache_items,
global_context,
)
-from dune.perftool.options import get_option, initialize_options
+from dune.perftool.options import (get_option,
+ initialize_options,
+ set_option_dependencies,
+ )
from dune.perftool.pdelab.driver import generate_driver
from dune.perftool.pdelab.localoperator import (generate_localoperator_basefile,
generate_localoperator_file,
@@ -113,6 +116,7 @@ def read_ufl(uflfile):
# This function is the entrypoint of the ufl2pdelab executable
def compile_form():
initialize_options()
+ set_option_dependencies()
formdatas, data = read_ufl(get_option("uflfile"))
with global_context(data=data, formdatas=formdatas):
diff --git a/python/dune/perftool/options.py b/python/dune/perftool/options.py
index 23733a7a..7afd80e1 100644
--- a/python/dune/perftool/options.py
+++ b/python/dune/perftool/options.py
@@ -42,6 +42,8 @@ class PerftoolOptionsArray(ImmutableRecord):
explicit_time_stepping = PerftoolOption(default=False, helpstr="use explicit time stepping")
exact_solution_expression = PerftoolOption(helpstr="name of the exact solution expression in the ufl file")
compare_l2errorsquared = PerftoolOption(helpstr="maximal allowed l2 error squared of difference between numerical solution and interpolation of exact solution (NOTE: requires --exact-solution-expression)")
+ l2error_tree_path = PerftoolOption(default=None, helpstr="Tree pathes that should be considered for l2 error calculation. Default None means we take all of them into account.")
+ interactive = PerftoolOption(default=False, helpstr="whether the optimization process should be guided interactively (also useful for debugging)")
print_transformations = PerftoolOption(default=False, helpstr="print out dot files after ufl tree transformations")
print_transformations_dir = PerftoolOption(default=".", helpstr="place where to put dot files (can be omitted)")
quadrature_order = PerftoolOption(_type=int, helpstr="Quadrature order used for all integrals.")
@@ -73,7 +75,8 @@ class PerftoolOptionsArray(ImmutableRecord):
precompute_quadrature_info = PerftoolOption(default=True, helpstr="compute quadrature points and weights in the constructor of the local operator")
blockstructured = PerftoolOption(default=False, helpstr="Use block structure")
number_of_blocks = PerftoolOption(default=1, helpstr="Number of sub blocks in one direction")
-
+ overlapping = PerftoolOption(default=False, helpstr="Use an overlapping solver and constraints. You still need to make sure to construct a grid with overlap! The parallel option will be set automatically.")
+ parallel = PerftoolOption(default=False, helpstr="Mark that this program should be run in parallel. If set to true the c++ code will check that there are more than 1 MPI-ranks involved and the error computation will use communication.")
# Until more sophisticated logic is needed, we keep the actual option data in this module
_options = PerftoolOptionsArray()
diff --git a/python/dune/perftool/pdelab/driver/__init__.py b/python/dune/perftool/pdelab/driver/__init__.py
index 04912744..14d0bd6c 100644
--- a/python/dune/perftool/pdelab/driver/__init__.py
+++ b/python/dune/perftool/pdelab/driver/__init__.py
@@ -251,10 +251,24 @@ def name_mpihelper():
return name
+@preamble
+def check_parallel_execution():
+ from dune.perftool.pdelab.driver.gridfunctionspace import name_leafview
+ gv = name_leafview()
+ return ["if ({}.comm().size()==1){{".format(gv),
+ ' std::cout << "This program should be run in parallel!" << std::endl;',
+ " return 1;",
+ "}"]
+
+
def generate_driver(formdatas, data):
# The driver module uses a global dictionary for storing necessary data
set_driver_data(formdatas, data)
+ # Add check to c++ file if this program should only be used in parallel mode
+ if get_option("parallel"):
+ check_parallel_execution()
+
# Entrypoint for driver generation
if get_option("opcounter") or get_option("time_opcounter"):
if get_option("time_opcounter"):
diff --git a/python/dune/perftool/pdelab/driver/constraints.py b/python/dune/perftool/pdelab/driver/constraints.py
index 9158fc6b..b6782295 100644
--- a/python/dune/perftool/pdelab/driver/constraints.py
+++ b/python/dune/perftool/pdelab/driver/constraints.py
@@ -29,9 +29,15 @@ def assemble_constraints(name):
element = get_trial_element()
gfs = name_trial_gfs()
is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- bctype_function = name_bctype_function(element, is_dirichlet)
- return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
- gfs,
+ from dune.perftool.pdelab.driver.interpolate import _do_interpolate
+ if _do_interpolate(is_dirichlet):
+ bctype_function = name_bctype_function(element, is_dirichlet)
+ return "Dune::PDELab::constraints({}, {}, {});".format(bctype_function,
+ gfs,
+ name,
+ )
+ else:
+ return "Dune::PDELab::constraints({}, {});".format(gfs,
name,
)
diff --git a/python/dune/perftool/pdelab/driver/error.py b/python/dune/perftool/pdelab/driver/error.py
index 6e7fb912..48bfe898 100644
--- a/python/dune/perftool/pdelab/driver/error.py
+++ b/python/dune/perftool/pdelab/driver/error.py
@@ -9,7 +9,8 @@ from dune.perftool.pdelab.driver import (get_formdata,
get_trial_element,
preprocess_leaf_data,
)
-from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_leafview,
+ name_trial_gfs,
name_trial_subgfs,
type_range,
)
@@ -104,14 +105,21 @@ def _accumulate_L2_squared(treepath):
strtp = ", ".join(str(t) for t in treepath)
+ gv = name_leafview()
+ sum_error_over_ranks = ""
+ if get_option("parallel"):
+ sum_error_over_ranks = " err = {}.comm().sum(err);".format(gv)
return ["{",
" // L2 error squared of difference between numerical",
" // solution and the interpolation of exact solution",
" // for treepath ({})".format(strtp),
" typename decltype({})::Traits::RangeType err(0.0);".format(dsa),
" Dune::PDELab::integrateGridFunction({}, err, 10);".format(dsa),
+ "{}".format(sum_error_over_ranks),
" {} += err;".format(accum_error),
- " std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+ " if ({}.comm().rank() == 0){{".format(gv),
+ " std::cout << \"L2 Error for treepath {}: \" << err << std::endl;".format(strtp),
+ " }"
"}",
]
@@ -139,8 +147,11 @@ def treepath_to_index(element, treepath, offset=0):
def accumulate_L2_squared():
element = get_trial_element()
if isinstance(element, MixedElement):
- for i in range(element.value_size()):
- _accumulate_L2_squared(get_treepath(element, i))
+ tree_pathes = list(map(int, get_option("l2error_tree_path").split(',')))
+ assert len(tree_pathes) == element.value_size()
+ for i, path in enumerate(tree_pathes):
+ if path:
+ _accumulate_L2_squared(get_treepath(element, i))
else:
_accumulate_L2_squared(())
@@ -160,12 +171,15 @@ def name_accumulated_L2_error():
@preamble
def compare_L2_squared():
accumulate_L2_squared()
+ gv = name_leafview()
accum_error = name_accumulated_L2_error()
fail = name_test_fail_variable()
return ["using std::abs;",
"using std::isnan;",
- "std::cout << \"\\nl2errorsquared: \" << {} << std::endl << std::endl;".format(accum_error),
+ "if ({}.comm().rank() == 0){{".format(gv),
+ " std::cout << \"\\nl2errorsquared: \" << {} << std::endl << std::endl;".format(accum_error),
+ "}",
"if (isnan({0}) or abs({0})>{1})".format(accum_error, get_option("compare_l2errorsquared")),
" {} = true;".format(fail)]
diff --git a/python/dune/perftool/pdelab/driver/gridfunctionspace.py b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
index f4392e36..c611f264 100644
--- a/python/dune/perftool/pdelab/driver/gridfunctionspace.py
+++ b/python/dune/perftool/pdelab/driver/gridfunctionspace.py
@@ -367,6 +367,18 @@ def type_orderingtag(leaf):
return "Dune::PDELab::EntityBlockedOrderingTag"
+@preamble
+def typedef_overlapping_dirichlet_constraintsassembler(name):
+ include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
+ return "using {} = Dune::PDELab::ConformingDirichletConstraints;".format(name)
+
+
+@preamble
+def typedef_p0parallel_constraintsassembler(name):
+ include_file("dune/pdelab/constraints/p0.hh", filetag="driver")
+ return "using {} = Dune::PDELab::P0ParallelConstraints;".format(name)
+
+
@preamble
def typedef_dirichlet_constraintsassembler(name):
include_file("dune/pdelab/constraints/conforming.hh", filetag="driver")
@@ -380,14 +392,21 @@ def typedef_no_constraintsassembler(name):
def type_constraintsassembler(is_dirichlet):
assert isinstance(is_dirichlet, bool)
- if is_dirichlet:
+ overlapping = get_option("overlapping")
+ if is_dirichlet and not overlapping:
name = "DirichletConstraintsAssember"
typedef_dirichlet_constraintsassembler(name)
- return name
+ elif is_dirichlet and overlapping:
+ name = "OverlappingConformingDirichletConstraints"
+ typedef_overlapping_dirichlet_constraintsassembler(name)
+ elif not is_dirichlet and overlapping:
+ name = "P0ParallelConstraints"
+ typedef_p0parallel_constraintsassembler(name)
else:
+ assert not is_dirichlet and not overlapping
name = "NoConstraintsAssembler"
typedef_no_constraintsassembler(name)
- return name
+ return name
def name_trial_subgfs(treepath):
diff --git a/python/dune/perftool/pdelab/driver/instationary.py b/python/dune/perftool/pdelab/driver/instationary.py
index 5c9bf3f6..ae4f3eb2 100644
--- a/python/dune/perftool/pdelab/driver/instationary.py
+++ b/python/dune/perftool/pdelab/driver/instationary.py
@@ -53,24 +53,37 @@ def time_loop():
params = name_parameters(formdata)
time = name_time()
element = get_trial_element()
- is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
- bctype = name_bctype_function(element, is_dirichlet)
- gfs = name_trial_gfs()
- cc = name_constraintscontainer()
vector_type = type_vector(formdata)
vector = name_vector(formdata)
interpolate_dirichlet_data(vector)
+ is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
+ from dune.perftool.pdelab.driver.interpolate import _do_interpolate
+ assemble_new_constraints = ""
+ dirichlet_constraints = _do_interpolate(is_dirichlet)
+ if dirichlet_constraints:
+ bctype = name_bctype_function(element, is_dirichlet)
+ gfs = name_trial_gfs()
+ cc = name_constraintscontainer()
+ assemble_new_constraints = (" // Assemble constraints for new time step\n"
+ " {}.setTime({}+dt);\n"
+ " Dune::PDELab::constraints({}, {}, {});\n"
+ "\n".format(params, time, bctype, gfs, cc)
+ )
+
# Choose between explicit and implicit time stepping
explicit = get_option('explicit_time_stepping')
if explicit:
osm = name_explicitonestepmethod()
apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
else:
+ osm = name_onestepmethod()
+ if dirichlet_constraints:
dirichlet = preprocess_leaf_data(element, "interpolate_expression")
boundary = name_boundary_function(element, dirichlet)
- osm = name_onestepmethod()
apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
+ else:
+ apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
# Setup visualization
visualize_initial_condition()
@@ -82,10 +95,7 @@ def time_loop():
"int step_number(0);"
"int output_every_nth = {}.get<int>(\"instat.output_every_nth\", 1);".format(ini),
"while (time<T-1e-8){",
- " // Assemble constraints for new time step",
- " {}.setTime({}+dt);".format(params, time),
- " Dune::PDELab::constraints({}, {}, {});".format(bctype, gfs, cc),
- "",
+ "{}".format(assemble_new_constraints),
" // Do time step",
" {} {}new({});".format(vector_type, vector, vector),
" {}".format(apply_call),
diff --git a/python/dune/perftool/pdelab/driver/solve.py b/python/dune/perftool/pdelab/driver/solve.py
index e877ece6..9b4cad91 100644
--- a/python/dune/perftool/pdelab/driver/solve.py
+++ b/python/dune/perftool/pdelab/driver/solve.py
@@ -7,7 +7,12 @@ from dune.perftool.pdelab.driver import (form_name_suffix,
is_linear,
name_initree,
)
-from dune.perftool.pdelab.driver.gridfunctionspace import name_trial_gfs
+from dune.perftool.pdelab.driver.gridfunctionspace import (name_trial_gfs,
+ type_trial_gfs,
+ )
+from dune.perftool.pdelab.driver.constraints import (type_constraintscontainer,
+ name_assembled_constraints,
+ )
from dune.perftool.pdelab.driver.gridoperator import (name_gridoperator,
type_gridoperator,
)
@@ -104,7 +109,12 @@ def define_vector(name, formdata):
@preamble
def typedef_linearsolver(name):
include_file("dune/pdelab/backend/istl.hh", filetag="driver")
- return "using {} = Dune::PDELab::ISTLBackend_SEQ_SuperLU;".format(name)
+ if get_option('overlapping'):
+ gfs = type_trial_gfs()
+ cc = type_constraintscontainer()
+ return "using {} = Dune::PDELab::ISTLBackend_OVLP_BCGS_ILU0<{},{}>;".format(name, gfs, cc)
+ else:
+ return "using {} = Dune::PDELab::ISTLBackend_SEQ_SuperLU;".format(name)
def type_linearsolver():
@@ -116,7 +126,12 @@ def type_linearsolver():
@preamble
def define_linearsolver(name):
lstype = type_linearsolver()
- return "{} {}(false);".format(lstype, name)
+ if get_option('overlapping'):
+ gfs = name_trial_gfs()
+ cc = name_assembled_constraints()
+ return "{} {}({}, {});".format(lstype, name, gfs, cc)
+ else:
+ return "{} {}(false);".format(lstype, name)
def name_linearsolver():
diff --git a/test/navier-stokes/howto/taylor-green.cc b/test/navier-stokes/howto/taylor-green.cc
index 218947d3..11150e51 100644
--- a/test/navier-stokes/howto/taylor-green.cc
+++ b/test/navier-stokes/howto/taylor-green.cc
@@ -58,7 +58,7 @@ void taylor_green (const GV& gv, const Dune::ParameterTree& configuration, std::
static const unsigned int dim = ES::dimension;
Dune::Timer timer;
- // Create grid function spaces
+ // Create finite element maps
const int velocity_degree = 2;
const int pressure_degree = 1;
typedef Dune::PDELab::QkDGLocalFiniteElementMap<DF, RF, velocity_degree, dim> FEM_V;
@@ -77,7 +77,6 @@ void taylor_green (const GV& gv, const Dune::ParameterTree& configuration, std::
using VectorBackend_V = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
using VectorBackend_P = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
using VectorBackend = Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::none>;
-
// typedef Dune::PDELab::EntityBlockedOrderingTag OrderingTag;
typedef Dune::PDELab::LexicographicOrderingTag OrderingTag_V;
@@ -254,21 +253,21 @@ void taylor_green (const GV& gv, const Dune::ParameterTree& configuration, std::
while (time < time_end - dt_min*0.5){
osm.apply(time,dt,xold,x);
- // Correct pressure after each step. Without this pressure
- // correction the velocity will be ok but the pressure will be
- // shifted by a constant.
- Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
- pressure_integral = gv.comm().sum(pressure_integral);
- std::cout << gv.comm().rank() << " palpo pressure_integral before: " << pressure_integral << std::endl;
-
- // Scale integral
- pressure_integral = pressure_integral/4;
- for (int i=pressure_index; i<gfs.size(); ++i){
- native(x)[i] -= pressure_integral;
- }
- Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
- pressure_integral = gv.comm().sum(pressure_integral);
- std::cout << "palpo pressure_integral after: " << pressure_integral << std::endl;
+ // // Correct pressure after each step. Without this pressure
+ // // correction the velocity will be ok but the pressure will be
+ // // shifted by a constant.
+ // Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
+ // pressure_integral = gv.comm().sum(pressure_integral);
+ // std::cout << gv.comm().rank() << " palpo pressure_integral before: " << pressure_integral << std::endl;
+
+ // // Scale integral
+ // pressure_integral = pressure_integral/4;
+ // for (int i=pressure_index; i<gfs.size(); ++i){
+ // native(x)[i] -= pressure_integral;
+ // }
+ // Dune::PDELab::integrateGridFunction(pdgf,pressure_integral,2);
+ // pressure_integral = gv.comm().sum(pressure_integral);
+ // std::cout << "palpo pressure_integral after: " << pressure_integral << std::endl;
xold = x;
time += dt;
diff --git a/test/navier-stokes/howto/taylor-green.ini b/test/navier-stokes/howto/taylor-green.ini
index 95250cc3..b7e94815 100644
--- a/test/navier-stokes/howto/taylor-green.ini
+++ b/test/navier-stokes/howto/taylor-green.ini
@@ -4,7 +4,7 @@ mu = 0.01
[parameters.dg]
epsilon = -1
-sigma = 10.0
+sigma = 6.0
beta = 1.0
[driver]
diff --git a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
index f916b0b9..bac1fc67 100644
--- a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
+++ b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.mini
@@ -14,10 +14,17 @@ extension = vtu
[formcompiler]
numerical_jacobian = 0, 1 | expand num
-compare_l2errorsquared = 1e-8
+compare_l2errorsquared = 5e-5
+# Only calculate error for the velocity part
+l2error_tree_path = 1, 1, 0
explicit_time_stepping = 0
grid_offset = 1
+overlapping = 1
[instat]
-T = 0.01
-dt = 0.0001
+T = 1e-2
+dt = 1e-3
+nth = 1
+
+# Disable numdiff tests
+{__exec_suffix} == numdiff | exclude
diff --git a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
index 39ab9c56..0797acc7 100644
--- a/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
+++ b/test/navier-stokes/navierstokes_2d_dg_quadrilateral.ufl
@@ -1,31 +1,35 @@
# Taylor-Green vortex
cell = quadrilateral
+degree = 2
+dim = 2
x = SpatialCoordinate(cell)
-# bctype = conditional(x[0] < 1. - 1e-8, 1, 0)
+time = variable(Constant(cell, count=2))
-P2 = VectorElement("DG", cell, 2)
-P1 = FiniteElement("DG", cell, 1)
+P2 = VectorElement("DG", cell, degree)
+P1 = FiniteElement("DG", cell, degree-1)
TH = P2 * P1
v, q = TestFunctions(TH)
u, p = TrialFunctions(TH)
-# ds = ds(subdomain_id=1, subdomain_data=bctype)
-
n = FacetNormal(cell)('+')
-eps = -1.0
-sigma = 1.0
-h_e = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
rho = 1.0
mu = 1.0/100.0
-t_start = 0.0
-g_v = as_vector((-exp(-2*pi*mu/rho*t_start)*cos(pi*x[0])*sin(pi*x[1]),
- exp(-2*pi*mu/rho*t_start)*sin(pi*x[0])*cos(pi*x[1])))
-g_p = -0.25*rho*exp(-4*pi*pi*mu/rho*t_start)*(cos(2*pi*x[0])+cos(2*pi*x[1]))
+g_v = as_vector((-exp(-2*pi*mu/rho*time)*cos(pi*x[0])*sin(pi*x[1]),
+ exp(-2*pi*mu/rho*time)*sin(pi*x[0])*cos(pi*x[1])))
+g_p = -0.25*rho*exp(-4*pi*pi*mu/rho*time)*(cos(2*pi*x[0])+cos(2*pi*x[1]))
+
+# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
+theta = -1.0
+
+# penalty factor
+alpha = 1.0
+h_int = Min(CellVolume(cell)('+'), CellVolume(cell)('-')) / FacetArea(cell)
+gamma_int = (alpha * degree * (degree + dim - 1)) / h_int
mass = rho*inner(u,v)*dx
@@ -34,16 +38,11 @@ r = mu * inner(grad(u), grad(v))*dx \
- q*div(u)*dx \
+ rho * inner(grad(u)*u,v)*dx \
+ mu * inner(avg(grad(u))*n, jump(v))*dS \
- + mu / h_e * sigma * inner(jump(u), jump(v))*dS \
- - mu * eps * inner(avg(grad(v))*n, jump(u))*dS \
+ + mu * gamma_int * inner(jump(u), jump(v))*dS \
+ - mu * theta * inner(avg(grad(v))*n, jump(u))*dS \
- avg(p)*inner(jump(v), n)*dS \
- avg(q)*inner(jump(u), n)*dS \
forms = [mass,r]
dirichlet_expression = g_v, g_p
-
-t_end = 1.0
-v_end = as_vector((-exp(-2*pi*mu/rho*t_end)*cos(pi*x[0])*sin(pi*x[1]),
- exp(-2*pi*mu/rho*t_end)*sin(pi*x[0])*cos(pi*x[1])))
-p_end = -0.25*rho*exp(-4*pi*pi*mu/rho*t_end)*(cos(2*pi*x[0])+cos(2*pi*x[1]))
-exact_solution = v_end, p_end
+exact_solution = g_v, g_p
--
GitLab