...
 
Commits (10)
#ifndef BLOCKSTRUCTURED_PRECONDITIONER_DECOMPOSITION_HH
#define BLOCKSTRUCTURED_PRECONDITIONER_DECOMPOSITION_HH
#include <algorithm>
#include <dune/pdelab/localoperator/flags.hh>
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/backend/istl.hh>
template<typename GFS, typename CC, typename DomainCount, int k>
template<typename GFS, typename CC, int k>
class SchurDecomposition{
using DF = double;
using RT = double;
class DomainCountLocalOperator: public Dune::PDELab::LocalOperatorDefaultFlags{
public:
static constexpr bool doAlphaVolume = true;
template<typename R, typename X, typename EG, typename LFSU, typename LFSV>
void alpha_volume(const EG &eg, const LFSU &lfsu, const X &x, const LFSV &lfsv, R &r) const {
auto& r_data = accessBaseContainer(r);
std::fill(begin(r_data), end(r_data), 1.);
}
};
using MB = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
using GOP = Dune::PDELab::GridOperator<GFS, GFS, DomainCount, MB, DF, RT, RT, CC, CC>;
using GOP = Dune::PDELab::GridOperator<GFS, GFS, DomainCountLocalOperator, MB, DF, RT, RT, CC, CC>;
using X = typename GOP::Domain;
using V = typename Dune::PDELab::Backend::Native<X>;
public:
explicit SchurDecomposition(const GFS& gfs, const CC& cc) {
using Dune::PDELab::Backend::native;
DomainCount lop(gfs, gfs, {});
DomainCountLocalOperator lop{};
GOP gop(gfs, cc, gfs, cc, lop, MB(1));
X zero(gfs, 0.), domainCount(gfs, 0.);
......
......@@ -75,11 +75,83 @@ def show_options():
print_options(form=True)
def process_global_options(opt):
""" Make sure that the options have been fully processed """
opt = expand_architecture_options(opt)
if opt.overlapping:
opt = opt.copy(parallel=True)
return opt
def process_form_options(opt, form):
if opt.sumfact:
opt = opt.copy(unroll_dimension_loops=True,
quadrature_mixins="sumfact",
basis_mixins="sumfact",
accumulation_mixins="sumfact",
)
if opt.blockstructured:
opt = opt.copy(accumulation_mixins="blockstructured",
quadrature_mixins="blockstructured",
basis_mixins="blockstructured"
)
if opt.control:
opt = opt.copy(accumulation_mixins="control")
if opt.numerical_jacobian:
opt = opt.copy(generate_jacobians=False, generate_jacobian_apply=False)
if opt.form is None:
opt = opt.copy(form=form)
if opt.classname is None:
opt = opt.copy(classname="{}Operator".format(form))
if opt.filename is None:
opt = opt.copy(filename="{}_{}_file.hh".format(get_option("target_name"), opt.classname))
if opt.block_preconditioner_pointdiagonal:
opt = opt.copy(generate_jacobians=False,
basis_mixins="sumfact_pointdiagonal",
accumulation_mixins="sumfact_pointdiagonal",
)
if opt.block_preconditioner_diagonal or opt.block_preconditioner_offdiagonal:
assert opt.numerical_jacobian is False
opt = opt.copy(generate_residuals=False,
generate_jacobians=True,
matrix_free=True,
)
if opt.matrix_free:
opt = opt.copy(generate_jacobian_apply=True)
return opt
def expand_architecture_options(opt):
if opt.architecture == "haswell":
return opt.copy(max_vector_width=256)
elif opt.architecture == "knl":
return opt.copy(max_vector_width=512)
elif opt.architecture == "skylake":
return opt.copy(max_vector_width=512)
else:
raise NotImplementedError("Architecture {} not known!".format(opt.architecture))
def initialize_options():
""" Initialize the options from the command line """
global _global_options
_global_options = update_options_from_commandline(_global_options)
_global_options = update_options_from_inifile(_global_options)
_global_options = process_global_options(_global_options)
for form in _form_options:
_form_options[form] = process_form_options(_form_options[form], form)
# Validate global options
scheme_global = _load_scheme()
......@@ -177,77 +249,6 @@ def update_options_from_inifile(opt):
return opt
@memoize
def process_global_options(opt):
""" Make sure that the options have been fully processed """
opt = expand_architecture_options(opt)
if opt.overlapping:
opt = opt.copy(parallel=True)
return opt
@memoize
def process_form_options(opt, form):
if opt.sumfact:
opt = opt.copy(unroll_dimension_loops=True,
quadrature_mixins="sumfact",
basis_mixins="sumfact",
accumulation_mixins="sumfact",
)
if opt.blockstructured:
opt = opt.copy(accumulation_mixins="blockstructured",
quadrature_mixins="blockstructured",
basis_mixins="blockstructured"
)
if opt.control:
opt = opt.copy(accumulation_mixins="control")
if opt.numerical_jacobian:
opt = opt.copy(generate_jacobians=False, generate_jacobian_apply=False)
if opt.form is None:
opt = opt.copy(form=form)
if opt.classname is None:
opt = opt.copy(classname="{}Operator".format(form))
if opt.filename is None:
opt = opt.copy(filename="{}_{}_file.hh".format(get_option("target_name"), opt.classname))
if opt.block_preconditioner_pointdiagonal:
opt = opt.copy(generate_jacobians=False,
basis_mixins="sumfact_pointdiagonal",
accumulation_mixins="sumfact_pointdiagonal",
)
if opt.block_preconditioner_diagonal or opt.block_preconditioner_offdiagonal:
assert opt.numerical_jacobian is False
opt = opt.copy(generate_residuals=False,
generate_jacobians=True,
matrix_free=True,
)
if opt.matrix_free:
opt = opt.copy(generate_jacobian_apply=True)
return opt
def expand_architecture_options(opt):
if opt.architecture == "haswell":
return opt.copy(max_vector_width=256)
elif opt.architecture == "knl":
return opt.copy(max_vector_width=512)
elif opt.architecture == "skylake":
return opt.copy(max_vector_width=512)
else:
raise NotImplementedError("Architecture {} not known!".format(opt.architecture))
def set_option(key, value):
"""Add the key value pair to the options.
......@@ -256,7 +257,7 @@ def set_option(key, value):
any other options.
"""
global _global_options
_global_options = process_global_options(_global_options).copy(**{key: value})
_global_options = _global_options.copy(**{key: value})
def set_form_option(key, value, form=None):
......@@ -269,8 +270,7 @@ def set_form_option(key, value, form=None):
def get_option(key):
processed_global_opts = process_global_options(_global_options)
return getattr(processed_global_opts, key)
return getattr(_global_options, key)
def get_form_option(key, form=None):
......@@ -279,8 +279,7 @@ def get_form_option(key, form=None):
form = get_global_context_value("form_identifier", 0)
if isinstance(form, int):
form = get_option("operators").split(",")[form].strip()
processed_form_opts = process_form_options(_form_options[form], form)
return getattr(processed_form_opts, key)
return getattr(_form_options[form], key)
@contextmanager
......
add_subdirectory(poisson)
add_subdirectory(nonlinear)
add_subdirectory(stokes)
\ No newline at end of file
add_subdirectory(stokes)
add_subdirectory(preconditioner)
\ No newline at end of file
......@@ -47,9 +47,3 @@ dune_add_formcompiler_system_test(UFLFILE poisson.ufl
BASENAME blockstructured_poisson_vec_tail
INIFILE poisson_vec_tail.mini
)
dune_add_formcompiler_system_test(UFLFILE poisson.ufl
BASENAME blockstructured_poisson_preconditioner
INIFILE poisson_preconditioner.mini
SOURCE preconditioner_driver.cc
)
\ No newline at end of file
dune_add_formcompiler_system_test(UFLFILE poisson.ufl
BASENAME blockstructured_preconditioner_poisson
INIFILE poisson.mini
SOURCE poisson_driver.cc
)
dune_add_formcompiler_system_test(UFLFILE linear_elasticity.ufl
BASENAME blockstructured_preconditioner_linear_elasticity
INIFILE linear_elasticity.mini
)
__name = blockstructured_preconditioner_linear_elasticity
cells = 40 2
extension = 1. 0.05
[wrapper.vtkcompare]
name = {__name}
extension = vtu
[formcompiler]
operators = r
[formcompiler.r]
blockstructured = 1
number_of_blocks = 8
matrix_free = 1
generate_jacobians = 0
geometry_mixins = blockstructured_equidistant
blockstructured_preconditioner = 1
[formcompiler.ufl_variants]
cell = quadrilateral
degree = 1
\ No newline at end of file
x = SpatialCoordinate(cell)
dim = as_domain(cell).geometric_dimension()
youngs_modulus = 200e9
poisson_coeff = 0.3
mu = 3 * youngs_modulus * (1 - 2 * poisson_coeff) / (2 * (1 + poisson_coeff))
lmd = 3 * youngs_modulus * poisson_coeff / (1 + poisson_coeff)
rho = 7800
g = 9.81
f = as_vector((0.,) * (dim - 1) + (-rho * g,))
bctype = conditional(abs(x[0]) < 1e-10, 1, 0)
V = VectorElement("CG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
def eps(u):
return 0.5 * (nabla_grad(u) + nabla_grad(u).T)
def sigma(u):
return lmd * nabla_div(u) * Identity(dim) + 2 * mu * eps(u)
r = (inner(sigma(u), eps(v)) - dot(f, v))*dx
interpolate_expression = (0.,) * dim
is_dirichlet = (bctype,) * dim
__name = blockstructured_poisson_preconditioner
__name = blockstructured_preconditioner_poisson
cells = 20 20
extension = 1. 1.
......
x = SpatialCoordinate(cell)
g = sum([c**2 for c in x])
f = sum([-2 for _ in x])
V = FiniteElement("CG", cell, degree)
u = TrialFunction(V)
v = TestFunction(V)
r = (inner(grad(u), grad(v)) - f*v)*dx
interpolate_expression = g
exact_solution = g
is_dirichlet = 1
......@@ -19,9 +19,8 @@
#include "dune/codegen/blockstructured/preconditioner/preconditioner.hh"
#include "dune/codegen/vtkpredicate.hh"
#include "blockstructured_poisson_preconditioner_rOperator_file.hh"
#include "blockstructured_preconditioner_poisson_rOperator_file.hh"
#include "rCoarseOperator_file.hh"
#include "domain_count_operator.hh"
#include "interpolation_local_operator.hh"
#include "restriction_local_operator.hh"
#include "local_decomposition_operator.hh"
......@@ -107,7 +106,7 @@ int main(int argc, char** argv){
auto func_0000 = Dune::PDELab::makeGridFunctionFromCallable(gv, lambda_0000);
Dune::PDELab::interpolate(func_0000, gfs, x_r);
SchurDecomposition<GFS, CC, DomainCountOperator<GFS, GFS>, NBLOCKS> decomp(gfs, cc);
SchurDecomposition<GFS, CC, NBLOCKS> decomp(gfs, cc);
using Decomposition = decltype(decomp);
RestrictionLocalOperator<GFS, GFS> r_lop(gfs, gfs, initree);
......