[feature][DGHelmholtzInv] Add possibility to pass update function for

preconditioner, i.e. set new linearization point.

dune/fem-dg/solver/dghelmholtzinverse.hh

@@ -16,6 +16,77 @@ namespace Dune
 {
   namespace Fem
   {
+    namespace detail
+    {
+      template <class DiscreteFunction>
+      class UpdatePreconditionerWrapper
+      {
+      public:
+        typedef DiscreteFunction DestinationType ;
+        typedef std::reference_wrapper< const DestinationType >  ConstReferenceType;
+        typedef std::function< void( ConstReferenceType& ) > UpdatePreconditionerFunctionType;
+
+        UpdatePreconditionerWrapper( const UpdatePreconditionerFunctionType& upPre )
+          : updatePrecond_( upPre )
+        {}
+
+        void update( const DestinationType& ubar ) const
+        {
+          ConstReferenceType u( ubar );
+          // call Python side update function
+          updatePrecond_( u );
+        }
+
+      protected:
+        const UpdatePreconditionerFunctionType& updatePrecond_;
+
+      };
+    }
+
+    template< class SpaceOperator >
+    class DGHelmholtzOperatorWithUpdate
+    : public DGHelmholtzOperator< SpaceOperator >
+    {
+      typedef DGHelmholtzOperatorWithUpdate< SpaceOperator > ThisType;
+      typedef DGHelmholtzOperator< SpaceOperator >           BaseType;
+
+    public:
+      typedef typename BaseType :: DomainFunctionType    DomainFunctionType;
+      typedef typename BaseType :: JacobianOperatorType  JacobianOperatorType;
+      typedef typename BaseType :: SpaceOperatorType     SpaceOperatorType ;
+
+      typedef detail::UpdatePreconditionerWrapper< DomainFunctionType > UpdatePreconditionerWrapperType;
+
+      using BaseType :: spaceOperator;
+      using BaseType :: lambda;
+
+      explicit DGHelmholtzOperatorWithUpdate ( SpaceOperatorType &spaceOp )
+        : BaseType( spaceOp )
+      {}
+
+      // overload jacobian to call update function
+      void jacobian ( const DomainFunctionType &u, JacobianOperatorType &jOp ) const
+      {
+        // update preconditioner to new linearization point
+        if( updatePreconditioner_ )
+          updatePreconditioner_->update( u );
+
+        spaceOperator().jacobian( u, jOp );
+        jOp.setLambda( lambda() );
+      }
+
+      void bind( const UpdatePreconditionerWrapperType& upPre )
+      {
+        updatePreconditioner_ = &upPre;
+      }
+
+      void unbind() { updatePreconditioner_ = nullptr; }
+
+    protected:
+      const UpdatePreconditionerWrapperType* updatePreconditioner_ = nullptr;
+    };
+
+
     /** \class DGHelmholtzInverseOperator
      *  \brief Operator implementing the solution of
      *
@@ -40,10 +111,13 @@ namespace Dune
       typedef std::reference_wrapper< DestinationType >        ReferenceType;
       typedef std::function< void( ConstReferenceType& , ReferenceType& ) > PreconditionerType;
 
+      typedef detail::UpdatePreconditionerWrapper< DestinationType > UpdatePreconditionerWrapperType;
+      typedef typename UpdatePreconditionerWrapperType :: UpdatePreconditionerFunctionType UpdatePreconditionerType;
+
     protected:
       typedef Dune::Fem::GmresInverseOperator< DestinationType >                LinearInverseOperatorType;
 
-      typedef Dune::Fem::DGHelmholtzOperator< SpaceOperatorType >               HelmholtzOperatorType;
+      typedef Dune::Fem::DGHelmholtzOperatorWithUpdate< SpaceOperatorType >     HelmholtzOperatorType;
       typedef Dune::Fem::NewtonInverseOperator< typename HelmholtzOperatorType::JacobianOperatorType,
                                                 LinearInverseOperatorType >     NonlinearInverseOperatorType;
 
@@ -75,6 +149,7 @@ namespace Dune
           pre_( uR, vR );
         }
       };
+      typedef PreconditionerWrapper PreconditionerWrapperType;
 
     public:
       DGHelmholtzInverseOperator( SpaceOperatorType& op,
@@ -122,6 +197,8 @@ namespace Dune
 
         // apply inv op
         (*this)( rhs, u );
+
+        helmholtzOp_.unbind();
         return SolverInfo( invOp_.converged(), invOp_.linearIterations(), invOp_.iterations() );
       }
 
@@ -132,14 +209,22 @@ namespace Dune
        * @f]
        *
        *  \param[in]     pre     preconditioner passed from Python side
+       *  \param[in]     upPre   update preconditioner function passed from Python side
        *  \param[in]     rhs     right hand side of F(w) = rhs
        *  \param[inout]  u       initial guess and returned solution
        *  \param[in]     lambda  lambda for Helmholtz operator
        */
       SolverInfo preconditionedSolve( const PreconditionerType& pre,
-                                      const DestinationType& rhs, DestinationType &u, const double lambda ) const
+                                      const UpdatePreconditionerType& upPre,
+                                      const DestinationType& rhs,
+                                      DestinationType &u,
+                                      const double lambda ) const
       {
-        PreconditionerWrapper p( pre );
+        PreconditionerWrapperType p( pre );
+        UpdatePreconditionerWrapperType up( upPre );
+
+        helmholtzOp_.bind( up );
+
         // bind operator and preconditioner
         invOp_.bind( helmholtzOp_, p );
         return solve( rhs, u, lambda );

python/dune/femdg/_operators.py

@@ -643,8 +643,12 @@ def dgHelmholtzInverseOperator( op, u = None, parameters = {} ):
                                  }''' )
     # add method solve, combining setLambda and __call__ for efficiency. Also,
     # here some solver diagnostics can be returned
-    preCondSolve = Method('preconditionedSolve', '''[]( DuneType &self, const typename DuneType::PreconditionerType& p, const typename DuneType::DestinationType &rhs, typename DuneType::DestinationType &u, const double lambda)
-                                 { auto info = self.preconditionedSolve(p, rhs, u, lambda);
+    preCondSolve = Method('preconditionedSolve', '''[]( DuneType &self,
+                                                        const typename DuneType::PreconditionerType& p,
+                                                        const typename DuneType::UpdatePreconditionerType& up,
+                                                        const typename DuneType::DestinationType &rhs,
+                                                        typename DuneType::DestinationType &u, const double lambda)
+                                 { auto info = self.preconditionedSolve(p, up, rhs, u, lambda);
                                    pybind11::dict ret;
                                    ret["converged"]  = pybind11::cast(info.converged);
                                    ret["iterations"] = pybind11::cast(info.nonlinearIterations);

python/dune/femdg/rk.py

@@ -155,9 +155,14 @@ class HelmholtzShuOsher:
         if self._invOp:
             # dummy preconditioner doing nothing
             #def pre(u,v):
+            #    # print("Pre: ", u.name, v.name )
             #    v.assign( u )
             #    return
-            #info = self._invOp.preconditionedSolve( pre, rhs, target, self.alpha )
+            # dummy updatePreconditioner doing nothing
+            #def updatePre( ubar ):
+            #    print("Update pre: ", ubar.name)
+            #    return
+            #info = self._invOp.preconditionedSolve( pre, updatePre, rhs, target, self.alpha )
 
             info = self._invOp.solve( rhs, target, self.alpha )
             self.counter = info["iterations"]