Skip to content
Snippets Groups Projects
Commit 757a6123 authored by René Heß's avatar René Heß
Browse files

Numerical matrix free solver for linear problems

parent 07c8b249
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/perftool/matrixfree.hh 0 → 100644
+ 353
− 0
View file @ 757a6123
#ifndef DUNE_PERFTOOL_MATRIXFREE_HH
#define DUNE_PERFTOOL_MATRIXFREE_HH
#include <iostream>
#include <dune/common/timer.hh>
#include <dune/common/parametertree.hh>
#include <dune/pdelab/backend/interface.hh>
#include <dune/pdelab/constraints/common/constraints.hh>
#include <dune/pdelab/backend/solver.hh>
namespace Dune{
namespace PDELab{
template <typename GO, typename V>
void solveMatrixFree(GO& go, V& x ){
typedef Dune::PDELab::OnTheFlyOperator<V,V,GO> ISTLOnTheFlyOperator;
ISTLOnTheFlyOperator opb(go);
Dune::Richardson<V,V> richardson(1.0);
Dune::CGSolver<V> solverb(opb,richardson,1E-10,5000,2);
Dune::InverseOperatorResult stat;
// evaluate residual w.r.t initial guess
using TrialGridFunctionSpace = typename GO::Traits::TrialGridFunctionSpace;
using W = Dune::PDELab::Backend::Vector<TrialGridFunctionSpace,typename V::ElementType>;
W r(go.testGridFunctionSpace(),0.0);
go.residual(x,r);
// solve the jacobian system
V z(go.trialGridFunctionSpace(),0.0);
solverb.apply(z,r,stat);
x -= z;
}
// template <template<class> class Solver, typename GO>
// class ISTLMatrixFreeSolverBackend_Base
// : public SequentialNorm, public LinearResultStorage
// {
// public:
// ISTLMatrixFreeSolverBackend_Base(GO go, unsigned maxiter_=5000, int verbose_=1)
// : _go(go), maxiter(maxiter_), verbose(verbose_)
// {}
// template <typename V, typename W>
// void apply(V z, W r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction) {
// using Backend::Native;
// using Backend::native;
// Dune::Richardson<V,V> prec(1.0);
// OnTheFlyOperator<V,V,GO> opa(_go);
// Solver<Native<V>> solver(opa, prec, reduction, maxiter, verbose);
// Dune::InverseOperatorResult stat;
// solver.apply(native(z), native(r), stat);
// }
// private:
// GO _go;
// unsigned maxiter;
// int verbose;
// };
// template <typename GO>
// class ISTLMatrixFreeBackend_SEQ_CG
// : public ISTLMatrixFreeSolverBackend_Base<Dune::CGSolver,GO>
// {
// public:
// ISTLMatrixFreeBackend_SEQ_CG (GO& go, unsigned maxiter_=5000, int verbose_=1)
// : ISTLMatrixFreeSolverBackend_Base<Dune::CGSolver,GO>(go, maxiter_, verbose_)
// {}
// };
// // Status information of linear problem solver
// template<class RFType>
// struct StationaryLinearMatrixFreeProblemSolverResult : LinearSolverResult<RFType>
// {
// RFType first_defect; // the first defect
// RFType defect; // the final defect
// double assembler_time; // Cumulative time for matrix assembly
// double linear_solver_time; // Cumulative time for linear solver
// int linear_solver_iterations; // Total number of linear iterations
// StationaryLinearMatrixFreeProblemSolverResult()
// : first_defect(0.0)
// , defect(0.0)
// , assembler_time(0.0)
// , linear_solver_time(0.0)
// , linear_solver_iterations(0)
// {}
// };
// template<typename GO, typename LS, typename V>
// class StationaryLinearMatrixFreeProblemSolver
// {
// typedef typename Dune::template FieldTraits<typename V::ElementType >::real_type Real;
// typedef typename GO::Traits::Jacobian M;
// typedef typename GO::Traits::TrialGridFunctionSpace TrialGridFunctionSpace;
// using W = Dune::PDELab::Backend::Vector<TrialGridFunctionSpace,typename V::ElementType>;
// typedef GO GridOperator;
// public:
// typedef StationaryLinearMatrixFreeProblemSolverResult<double> Result;
// StationaryLinearMatrixFreeProblemSolver(GO& go, LS& ls, V& x, Real reduction, Real min_defect = 1e-99, int verbose=1)
// : _go(go)
// , _ls(ls)
// , _x(&x)
// , _reduction(reduction)
// , _min_defect(min_defect)
// , _hanging_node_modifications(false)
// , _keep_matrix(true)
// , _verbose(verbose)
// {}
// StationaryLinearMatrixFreeProblemSolver (const GO& go, LS& ls, Real reduction, Real min_defect = 1e-99, int verbose=1)
// : _go(go)
// , _ls(ls)
// , _x()
// , _reduction(reduction)
// , _min_defect(min_defect)
// , _hanging_node_modifications(false)
// , _keep_matrix(true)
// , _verbose(verbose)
// {}
// //! Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterTree.
// /**
// * This constructor reads the parameter controlling its operation from a passed-in ParameterTree
// * instead of requiring the user to specify all of them as individual constructor parameters.
// * Currently the following parameters are read:
// *
// * Name | Default Value | Explanation
// * -------------------------- | ------------- | -----------
// * reduction | | Required relative defect reduction
// * min_defect | 1e-99 | minimum absolute defect at which to stop
// * hanging_node_modifications | false | perform required transformations for hanging nodes
// * keep_matrix | true | keep matrix between calls to apply() (but reassemble values every time)
// * verbosity | 1 | control amount of debug output
// *
// * Apart from reduction, all parameters have a default value and are optional.
// * The actual reduction for a call to apply() is calculated as r = max(reduction,min_defect/start_defect),
// * where start defect is the norm of the residual of x.
// */
// StationaryLinearMatrixFreeProblemSolver(const GO& go, LS& ls, V& x, const ParameterTree& params)
// : _go(go)
// , _ls(ls)
// , _x(&x)
// , _reduction(params.get<Real>("reduction"))
// , _min_defect(params.get<Real>("min_defect",1e-99))
// , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
// , _keep_matrix(params.get<bool>("keep_matrix",true))
// , _verbose(params.get<int>("verbosity",1))
// {}
// //! Construct a StationaryLinearProblemSolver for the given objects and read parameters from a ParameterTree.
// /**
// * This constructor reads the parameter controlling its operation from a passed-in ParameterTree
// * instead of requiring the user to specify all of them as individual constructor parameters.
// * Currently the following parameters are read:
// *
// * Name | Default Value | Explanation
// * -------------------------- | ------------- | -----------
// * reduction | | Required relative defect reduction
// * min_defect | 1e-99 | minimum absolute defect at which to stop
// * hanging_node_modifications | false | perform required transformations for hanging nodes
// * keep_matrix | true | keep matrix between calls to apply() (but reassemble values every time)
// * verbosity | 1 | control amount of debug output
// *
// * Apart from reduction, all parameters have a default value and are optional.
// * The actual reduction for a call to apply() is calculated as r = max(reduction,min_defect/start_defect),
// * where start defect is the norm of the residual of x.
// */
// StationaryLinearMatrixFreeProblemSolver(const GO& go, LS& ls, const ParameterTree& params)
// : _go(go)
// , _ls(ls)
// , _x()
// , _reduction(params.get<Real>("reduction"))
// , _min_defect(params.get<Real>("min_defect",1e-99))
// , _hanging_node_modifications(params.get<bool>("hanging_node_modifications",false))
// , _keep_matrix(params.get<bool>("keep_matrix",true))
// , _verbose(params.get<int>("verbosity",1))
// {}
// // //! Set whether the solver should apply the necessary transformations for calculations on hanging nodes.
// // void setHangingNodeModifications(bool b)
// // {
// // _hanging_node_modifications=b;
// // }
// // //! Return whether the solver performs the necessary transformations for calculations on hanging nodes.
// // bool hangingNodeModifications() const
// // {
// // return _hanging_node_modifications;
// // }
// // //! Set whether the jacobian matrix should be kept across calls to apply().
// // void setKeepMatrix(bool b)
// // {
// // _keep_matrix = b;
// // }
// // //! Return whether the jacobian matrix is kept across calls to apply().
// // bool keepMatrix() const
// // {
// // return _keep_matrix;
// // }
// const Result& result() const
// {
// return _res;
// }
// void apply(V& x, bool reuse_matrix = false) {
// _x = &x;
// apply(reuse_matrix);
// }
// void apply (bool reuse_matrix = false)
// {
// Dune::Timer watch;
// double timing,assembler_time=0;
// // // assemble matrix; optional: assemble only on demand!
// // watch.reset();
// // if (!_jacobian)
// // {
// // _jacobian = std::make_shared<M>(_go);
// // timing = watch.elapsed();
// // if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
// // std::cout << "=== matrix setup (max) " << timing << " s" << std::endl;
// // watch.reset();
// // assembler_time += timing;
// // }
// // else if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
// // std::cout << "=== matrix setup skipped (matrix already allocated)" << std::endl;
// // if (_hanging_node_modifications)
// // {
// // Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
// // _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
// // }
// // if (!reuse_matrix)
// // {
// // (*_jacobian) = Real(0.0);
// // _go.jacobian(*_x,*_jacobian);
// // }
// // timing = watch.elapsed();
// // // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
// // if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
// // {
// // if (reuse_matrix)
// // std::cout << "=== matrix assembly SKIPPED" << std::endl;
// // else
// // std::cout << "=== matrix assembly (max) " << timing << " s" << std::endl;
// // }
// // assembler_time += timing;
// // assemble residual
// watch.reset();
// W r(_go.testGridFunctionSpace(),0.0);
// _go.residual(*_x,r); // residual is additive
// timing = watch.elapsed();
// // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
// if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
// std::cout << "=== residual assembly (max) " << timing << " s" << std::endl;
// assembler_time += timing;
// _res.assembler_time = assembler_time;
// auto defect = _ls.norm(r);
// // compute correction
// watch.reset();
// V z(_go.trialGridFunctionSpace(),0.0);
// auto red = std::max(_reduction,_min_defect/defect);
// if (_go.trialGridFunctionSpace().gridView().comm().rank()==0)
// std::cout << "=== solving (reduction: " << red << ") " << std::endl;
// // NOTE: Different than in StationaryLinearProblemSolver!
// _ls.apply(z, r, red); // solver makes right hand side consistent
// _linear_solver_result = _ls.result();
// timing = watch.elapsed();
// // timing = gos.trialGridFunctionSpace().gridView().comm().max(timing);
// if (_go.trialGridFunctionSpace().gridView().comm().rank()==0 && _verbose>=1)
// std::cout << timing << " s" << std::endl;
// _res.linear_solver_time = timing;
// _res.converged = _linear_solver_result.converged;
// _res.iterations = _linear_solver_result.iterations;
// _res.elapsed = _linear_solver_result.elapsed;
// _res.reduction = _linear_solver_result.reduction;
// _res.conv_rate = _linear_solver_result.conv_rate;
// _res.first_defect = static_cast<double>(defect);
// _res.defect = static_cast<double>(defect)*_linear_solver_result.reduction;
// _res.linear_solver_iterations = _linear_solver_result.iterations;
// // // and update
// // if (_hanging_node_modifications)
// // Dune::PDELab::set_shifted_dofs(_go.localAssembler().trialConstraints(),0.0,*_x); // set hanging node DOFs to zero
// // *_x -= z;
// // if (_hanging_node_modifications)
// // _go.localAssembler().backtransform(*_x); // interpolate hanging nodes adjacent to Dirichlet nodes
// // if (!_keep_matrix)
// // _jacobian.reset();
// }
// // //! Discard the stored Jacobian matrix.
// // void discardMatrix()
// // {
// // if(_jacobian)
// // _jacobian.reset();
// // }
// const Dune::PDELab::LinearSolverResult<double>& ls_result() const{
// return _linear_solver_result;
// }
// Real reduction() const
// {
// return _reduction;
// }
// void setReduction(Real reduction)
// {
// _reduction = reduction;
// }
// private:
// const GO& _go;
// LS& _ls;
// V* _x;
// // shared_ptr<M> _jacobian;
// Real _reduction;
// Real _min_defect;
// Dune::PDELab::LinearSolverResult<double> _linear_solver_result;
// Result _res;
// bool _hanging_node_modifications;
// bool _keep_matrix;
// int _verbose;
// };
}
}
#endif
python/dune/perftool/options.py
+ 1
− 0
View file @ 757a6123
......@@ -13,6 +13,7 @@ def get_form_compiler_arguments():
parser.add_argument("--operator-file", type=str, help="The filename for the generated local operator header")
parser.add_argument('--version', action='version', version='%(prog)s 0.1')
parser.add_argument("--numerical-jacobian", action="store_true", help="use numerical jacobians (only makes sense, if uflpdelab for some reason fails to generate analytic jacobians)")
parser.add_argument("--matrix-free", action="store_true", help="use iterative solver with matrix free jacobian application")
parser.add_argument(
"--exact-solution-expression",
type=str,
......
python/dune/perftool/pdelab/driver.py
+ 8
− 2
View file @ 757a6123
......@@ -903,8 +903,14 @@ def name_stationarynonlinearproblemsolver():
def dune_solve():
# Test if form is linear in ansatzfunction
if is_linear(_form):
slp = name_stationarylinearproblemsolver()
return "{}.apply();".format(slp)
if not get_option("matrix_free"):
slp = name_stationarylinearproblemsolver()
return "{}.apply();".format(slp)
else:
go = name_gridoperator()
x = name_vector()
include_file("dune/perftool/matrixfree.hh", filetag="drive")
return "solveMatrixFree({},{});".format(go,x)
else:
snp = name_stationarynonlinearproblemsolver()
return "{}.apply();".format(snp)
......
python/dune/perftool/pdelab/localoperator.py
+ 14
− 0
View file @ 757a6123
......@@ -319,16 +319,30 @@ def generate_localoperator_kernels(formdata, data):
_, loptype = class_type_from_cache("operator")
which = ufl_measure_to_pdelab_measure(measure)
# Numerical jacobian methods
base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
# Add the initializer list for that base class
ini = name_initree_member()
ini_constructor = name_initree_constructor_param()
# Numerical jacobian methods
initializer_list("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
classtag="operator",
)
if get_option("matrix_free"):
# Numeical jacobian apply base class
base_class("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype), classtag="operator")
# Numerical jacobian apply initializer list
initializer_list("Dune::PDELab::NumericalJacobianApply{}<{}>".format(which, loptype),
["{}.get<double>(\"numerical_epsilon.{}\", 1e-9)".format(ini_constructor, ini, which.lower())],
classtag="operator",
)
operator_kernels[(measure, 'residual')] = kernel
# Generate the necessary jacobian methods
......
test/poisson/CMakeLists.txt
+ 12
− 0
View file @ 757a6123
......@@ -88,6 +88,18 @@ dune_add_system_test(TARGET poisson_dg_neumann_symdiff
INIFILE poisson_dg_neumann_symdiff.mini
)
# 9. Poisson Test Case: source f, bc g + numerical differentiation
add_generated_executable(UFLFILE poisson.ufl
TARGET poisson_numdiff_matrix_free
FORM_COMPILER_ARGS --numerical-jacobian --matrix-free
)
dune_add_system_test(TARGET poisson_numdiff_matrix_free
INIFILE poisson_numdiff_matrix_free.mini
SCRIPT dune_vtkcompare.py
)
# Add the reference code with the PDELab localoperator that produced
# the reference vtk file
add_executable(poisson_dg_ref reference_main.cc)
......
test/poisson/poisson_numdiff_matrix_free.mini 0 → 100644
+ 11
− 0
View file @ 757a6123
__name = poisson_numdiff_matrix_free
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 32 32
elementType = simplical
[wrapper.vtkcompare]
name = {__name}
reference = poisson_ref
extension = vtu
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Datenschutz und Impressum