Commit 74a94dd2 authored by Marcel Koch's avatar Marcel Koch

adapt to pdelab changes

parent 86ec58f2
Pipeline #25362 failed with stage
in 24 minutes and 4 seconds
......@@ -11,8 +11,8 @@
namespace Dune{
namespace PDELab{
template <typename GO, typename V, typename Preconditioner = Dune::Richardson<V,V>>
void solveMatrixFree(GO& go, V& x, Preconditioner& pre = Preconditioner{1.0}){
template <typename GO, typename V, typename Preconditioner>
void solveMatrixFree(GO& go, V& x, Preconditioner& pre){
using ISTLOnTheFlyOperator = Dune::PDELab::OnTheFlyOperator<V,V,GO>;
ISTLOnTheFlyOperator opb(go);
Dune::BiCGSTABSolver<V> solverb(opb,pre,1E-10,5000,2);
......@@ -28,16 +28,26 @@ namespace Dune{
x -= z;
}
template <typename GO, typename V>
void solveMatrixFree(GO& go, V& x){
using LinearSolver = ISTLBackend_SEQ_MatrixFree_BCGS_Richardson<GO>;
LinearSolver ls(go, 5000, 2);
// 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);
ls.apply(z,r,1e-10);
x -= z;
}
template <typename GO, typename V>
void solveNonlinearMatrixFree(GO& go, V& x ){
// Matrix free operator
using ISTLOnTheFlyOperator = Dune::PDELab::NonlinearOnTheFlyOperator<V,V,GO>;
ISTLOnTheFlyOperator op(go);
// Matrix free linear solver
using LinearSolver = ISTLBackend_SEQ_MatrixFree_Richardson<ISTLOnTheFlyOperator, Dune::BiCGSTABSolver>;
LinearSolver ls(op);
using LinearSolver = ISTLBackend_SEQ_MatrixFree_BCGS_Richardson<GO>;
LinearSolver ls(go);
// Solve nonlinear system with matrix free newton
using SNP = MatrixFreeNewton<GO, LinearSolver, V>;
......
......@@ -39,6 +39,7 @@ from dune.codegen.cgen.clazz import (AccessModifier,
BaseClass,
ClassMember,
)
from dune.codegen.pdelab.driver import is_linear
from dune.codegen.loopy.target import type_floatingpoint
from dune.codegen.ufl.modified_terminals import Restriction
from frozendict import frozendict
......@@ -157,6 +158,15 @@ def enum_skeleton_twosided():
return "enum { doSkeletonTwoSided = true };"
@class_member(classtag="operator")
def enum_is_linear():
linear = is_linear()
if linear:
return "enum { isLinear = true };"
else:
return "enum { isLinear = false };"
def name_initree_member():
define_initree("_iniParams")
return "_iniParams"
......@@ -919,6 +929,9 @@ def local_operator_default_settings(operator, form):
enum_pattern()
pattern_baseclass()
# For non-linear local operators we need to set the isLinear flag to false
enum_is_linear()
if get_form_option("block_preconditioner_diagonal") or get_form_option("block_preconditioner_pointdiagonal"):
enum_skeleton_twosided()
......
Markdown is supported
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