mass conserving problemstat not using multi-mesh added

extensions/ProblemStatMassConserve2.h

+/******************************************************************************
+ *
+ * Mass-conserving ProblemStat for AMDiS
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis/trunk/extensions
+ *
+ * Author: Simon Praetorius
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+/** \file ProblemStatMassConserve2.h */
+
+
+#ifndef PROBLEM_STAT_MASS_CONSERVE
+#define PROBLEM_STAT_MASS_CONSERVE
+
+#include "ExtendedProblemStat.h"
+#include "Refinement.h"
+
+/// Implementation of ProblemStat to allow for the conservation of mass of one
+/// solution component, i.e. project the old-solution to the new mesh so that
+/// the scalar-product is conserved: \f$ (u^*, v) = (u^\text{old}, v) \f$, where \f$ u^*,v \f$ live on the
+/// new mesh and \f$ u^\text{old} \f$ lives on the old mesh.
+struct ProblemStatMassConserve2 : public ExtendedProblemStat, public StandardRefineOperation
+{
+  /// constructor: create a temporary problem that is initialized with 
+  /// similar parameters as the regular one
+  ProblemStatMassConserve2(std::string name_)
+  : ExtendedProblemStat(name_),
+    prob2(nullptr),
+    comp(0)
+  {    
+    Parameters::get(name_ + "->mass conservation component", comp);
+    
+    // set parameters for the temporary problem
+    std::string name2 = name_ + "_tmp";
+    int zero = 0;
+    int one = 1;
+    int two = 2;
+    double tol = 1.e-8;
+    int degree = 1;
+    int dim = 1;
+    std::string meshName = "";
+    
+    Parameters::get(name_ + "->polynomial degree[" + boost::lexical_cast<std::string>(comp) + "]", degree);
+    Parameters::get(name_ + "->dim", dim);
+    Parameters::get(name_ + "->mesh", meshName);
+    
+    Parameters::set(name2 + "->components", one);
+    Parameters::set(name2 + "->polynomial degree[0]", degree);
+    Parameters::set(name2 + "->dim", dim);
+    Parameters::set(name2 + "->mesh", meshName);
+#if defined HAVE_UMFPACK || defined HAVE_PARALLEL_PETSC || defined HAVE_PETSC
+    std::string solverName = "direct";
+#else
+    std::string solverName = "cg";
+#endif
+    Parameters::set(name2 + "->solver", solverName);
+    Parameters::set(name2 + "->tolerance", tol);
+    Parameters::set(name2 + "->info", zero);
+    
+    // create a new temporary problem
+    prob2 = new ProblemStat(name2);
+  }
+  
+  
+  ~ProblemStatMassConserve2()
+  {
+    if (prob2 != nullptr) {
+      delete prob2;
+      prob2 = nullptr;
+    }
+  }
+  
+  
+  /// make shure that temporary problem has an equal mesh as regular problem
+  /// and initialize this mesh only once
+  void initialize(Flag initFlag,
+		  ProblemStatSeq *adoptProblem = NULL,
+		  Flag adoptFlag = INIT_NOTHING) override
+  {
+    ExtendedProblemStat::initialize(initFlag - INIT_MESH, adoptProblem, adoptFlag); // create meshes for prob
+    
+    std::vector<Mesh*> meshes2;
+    meshes2.push_back(getMesh(comp));
+    prob2->setComponentMesh(0, getMesh(comp));
+    prob2->setMeshes(meshes2);    
+    prob2->initialize(INIT_ALL - INIT_MESH);
+    
+    for (size_t i = 0; i < meshes2.size(); i++) {
+      int globalRefinements = 0;
+
+      // If AMDiS is compiled for parallel computations, the global refinements are
+      // ignored here. Later, each rank will add the global refinements to its
+      // private mesh.
+#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
+      Parameters::get(meshes2[i]->getName() + "->global refinements",
+		      globalRefinements);
+#endif
+
+      bool initMesh = true;
+
+      // Initialize the meshes if there is no serialization file.
+      if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) {
+	meshes2[i]->initialize();
+	prob2->getRefinementManager()->globalRefine(meshes2[i], globalRefinements);
+      }
+    }
+    
+    prob2->getRhsVector(comp)->setCoarsenOperation(COARSE_RESTRICT);
+    
+    fillOperators();
+  }
+  
+  /// implementation of StandardRefineOperation::beforeCoarsen
+  void beforeCoarsen(AdaptInfo* adaptInfo, Flag markFlag) override
+  {
+    if (markFlag.isSet(MESH_COARSENED))
+      prob2->buildAfterCoarsen(adaptInfo, markFlag, false, true);
+  }
+  
+  /// in the first iteration copy initial-solution to temporary problem
+  /// and prepare the temporary mesh
+  void buildBeforeCoarsen(AdaptInfo* adaptInfo, Flag markFlag) override
+  {   
+    ExtendedProblemStat::buildBeforeCoarsen(adaptInfo, markFlag); 
+    beforeRefine(adaptInfo, markFlag);
+  }
+  
+  
+  /// implementation of StandardRefineOperation::afterCoarsen
+  void afterCoarsen(AdaptInfo* adaptInfo, Flag markFlag) override
+  {
+    if (markFlag.isSet(MESH_COARSENED)) {
+      prob2->buildAfterCoarsen(adaptInfo, markFlag, true, false);
+      prob2->solve(adaptInfo);
+      *getSolution(comp) = *prob2->getSolution(0);
+    }
+  }
+  
+  /// solve temporary problem to project solution to new mesh and then
+  /// copy temporary solution back to regular solution
+  void buildAfterCoarsen(AdaptInfo* adaptInfo, 
+			 Flag markFlag, 
+			 bool asmM, 
+			 bool asmV) override
+  {  
+    afterCoarsen(adaptInfo, markFlag);
+    
+    // call original buildAfterCoarsen
+    ExtendedProblemStat::buildAfterCoarsen(adaptInfo, markFlag, asmM, asmV); 
+  }
+  
+  /// add mass-matrix and transfer-operator to temporary problem, to
+  /// project solution to new mesh
+  void fillOperators()
+  {
+    // conversion propblem 2
+    const FiniteElemSpace* feSpace2 = prob2->getFeSpace(0);
+
+    Operator *opMnew0 = new Operator(feSpace2, feSpace2);
+    opMnew0->addTerm(new Simple_ZOT);
+    Operator *opMold01 = new Operator(feSpace2, feSpace2);
+    opMold01->addTerm(new VecAtQP_ZOT(getSolution(comp)));
+    
+    // assemble only one block of the system
+    prob2->addMatrixOperator(*opMnew0,0,0);
+    prob2->addVectorOperator(*opMold01, 0); 
+  }
+  
+private:  
+  /// temporary problem for the projection of the old-solution to the new mesh
+  ProblemStat* prob2;
+  
+  /// component the should be copnserved
+  int comp;
+};
+
+#endif // PROBLEM_STAT_MASS_CONSERVE
\ No newline at end of file

extensions/Refinement.h

@@ -86,6 +86,14 @@ protected:
 };
 
 
+struct StandardRefineOperation
+{
+  virtual void beforeRefine(AdaptInfo* adaptInfo, Flag markFlag) {}
+  virtual void beforeCoarsen(AdaptInfo* adaptInfo, Flag markFlag) {}
+  virtual void afterCoarsen(AdaptInfo* adaptInfo, Flag markFlag) {}  
+};
+
+
 /** \brief
  * Base class for Refinement structure to perform local anisotropic refinement
  */
@@ -94,9 +102,12 @@ class RefinementLevel
 {
 public:
 
-  RefinementLevel(const FiniteElemSpace *feSpace_, MeshRefinementFunction<T,T2>* refineFct_) :
+  RefinementLevel(const FiniteElemSpace *feSpace_, 
+		  MeshRefinementFunction<T,T2>* refineFct_) :
     feSpace(feSpace_), 
     refineFct(refineFct_),
+    adaptInfo(nullptr),
+    refineOperation(nullptr),
     numRefinements0(15),
     globalRefined(false)
   {	
@@ -121,11 +132,25 @@ public:
     }
 
     numRefinements = numRefinements0;
+    
+    refineOperation = new StandardRefineOperation;
   }
   
   virtual ~RefinementLevel() {
     delete coarseningManager;
     delete refinementManager;
+    if (refineOperation)
+      delete refineOperation;
+  }
+  
+  void setRefineOperation(AdaptInfo* adaptInfo_,
+			  StandardRefineOperation* refineOperation_)
+  {
+    if (refineOperation)
+      delete refineOperation;
+    
+    adaptInfo = adaptInfo_;
+    refineOperation = refineOperation_;
   }
 
   void refine(bool onlyRefine= false) 
@@ -147,7 +172,7 @@ public:
     int i = 0;
     while (meshChanged && i < numRefinements) {
       markElements(markFlag);
-      meshChanged = refineMesh(markFlag, onlyRefine);
+      meshChanged = refineMesh(adaptInfo, markFlag, onlyRefine);
       
       calcMeshSizes(minH, maxH, minLevel, maxLevel, false); 
       int nr = mesh->getNumberOfVertices();
@@ -258,10 +283,13 @@ public:
   bool refineMesh(Flag markFlag, bool onlyRefine) 
   {
     int oldSize = mesh->getNumberOfVertices();
+    refineOperation->beforeRefine(adaptInfo, markFlag);
     if (markFlag.isSet(1))
       refinementManager->refineMesh(mesh);
+    refineOperation->beforeCoarsen(adaptInfo, markFlag);
     if (markFlag.isSet(2) && !onlyRefine)
       coarseningManager->coarsenMesh(mesh);
+    refineOperation->afterCoarsen(adaptInfo, markFlag);
     if (markFlag.isSet(1) || markFlag.isSet(2)) {
       int newSize = mesh->getNumberOfVertices();
       if (oldSize != newSize) 
@@ -283,9 +311,13 @@ protected:
   MeshRefinementFunction<T,T2>* refineFct;
   RefinementManager* refinementManager;
   CoarseningManager* coarseningManager;
+  
+  AdaptInfo* adaptInfo;
+  StandardRefineOperation* refineOperation;
 
   int numRefinements;
   int numRefinements0;
+  
   bool globalRefined;
 };
 

extensions/base_problems/BaseProblem.hh

@@ -268,7 +268,8 @@ Flag BaseProblem<ProblemType, safeguard>::oneIteration(AdaptInfo *adaptInfo,
   if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
     flag = prob->refineMesh(adaptInfo);
 
-  prob->buildBeforeCoarsen(adaptInfo, markFlag);
+  if (toDo.isSet(BUILD))
+    prob->buildBeforeCoarsen(adaptInfo, markFlag);
 
   // coarsen
   if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))

extensions/demo/other/init/cahnHilliard_navierStokes.dat.2d

@@ -39,11 +39,15 @@ ns->non-linear term: 2
 %   3... (1)+(2)-(0)
 
 
+
+ch->space->mass conservation component: 0
+
 % =========== OUTPUT ==============================================
-ch->space->output->filename: ./output/ch_
+prefix: _restrict2
+ch->space->output->filename: ./output/ch${prefix}_
 
-ns->velocity->output->filename: ./output/velocity_
-ns->space->output[2]->filename: ./output/pressure_
+ns->velocity->output->filename: ./output/velocity${prefix}_
+ns->space->output[2]->filename: ./output/pressure${prefix}_
 
 % ==================== TIMESTEPS ===============================
 adapt->timestep:      1.e-2

extensions/demo/other/init/elliptBase.dat.2d

+dimension of world:             2
+
+elliptMesh->macro file name:    ./macro/stand.macro.2d
+elliptMesh->global refinements: 5
+
+ellipt->space->mesh:                   elliptMesh
+ellipt->space->dim:                    2
+ellipt->space->components:             1
+ellipt->space->polynomial degree[0]:   1
+ 
+ellipt->space->solver:                 cg
+ellipt->space->solver->ell:        1
+ellipt->space->solver->max iteration:  1000
+ellipt->space->solver->tolerance:      1.e-8
+ellipt->space->solver->info:           10
+ellipt->space->solver->left precon:    ilu
+ellipt->space->solver->right precon:   no
+
+ellipt->space->estimator[0]:              residual
+ellipt->space->estimator[0]->error norm:  H1_NORM   % 1: H1_NORM, 2: L2_NORM
+ellipt->space->estimator[0]->C0:          0.1 % constant of element residual
+ellipt->space->estimator[0]->C1:          0.1 % constant of jump residual
+
+ellipt->space->marker[0]->strategy:       2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+ellipt->space->marker[0]->MSGamma:        0.5
+
+ellipt->adapt[0]->tolerance:          1e-6
+ellipt->adapt[0]->refine bisections:  2
+
+ellipt->adapt->max iteration:         10
+
+ellipt->space->output->filename:       output/elliptBase.2d
+ellipt->space->output->ParaView format: 1
+ellipt->space->output->ARH format: 1
+ellipt->space->output->AMDiS format: 1

extensions/demo/other/init/elliptImplicit.dat.2d

+dimension of world:             2
+
+elliptMesh->macro file name:    ./macro/macro.square.2d
+elliptMesh->global refinements: 0
+
+mesh->refinement->initial level: 10
+mesh->refinement->level on interface: 17
+mesh->refinement->level in inner domain: 10
+mesh->refinement->level in outer domain: 0
+mesh->refinement->interface width: 0.01
+
+ellipt->mesh:                   elliptMesh
+ellipt->dim:                    2
+ellipt->components:             1
+ellipt->polynomial degree[0]:   2
+
+ellipt->solver:                 cg
+ellipt->solver->ell:        1
+ellipt->solver->max iteration:  1000
+ellipt->solver->tolerance:      1.e-8
+ellipt->solver->info:           10
+ellipt->solver->left precon:    ilu
+ellipt->solver->right precon:   no
+
+ellipt->estimator[0]:              residual
+ellipt->estimator[0]->error norm:  H1_NORM   % 1: H1_NORM, 2: L2_NORM
+ellipt->estimator[0]->C0:          0.1 % constant of element residual
+ellipt->estimator[0]->C1:          0.1 % constant of jump residual
+
+ellipt->marker[0]->strategy:       2 % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+ellipt->marker[0]->MSGamma:        0.5
+
+ellipt->adapt[0]->tolerance:          1e-6
+ellipt->adapt[0]->refine bisections:  2
+
+ellipt->adapt->max iteration:         10
+
+ellipt->output->filename:       output/elliptImplicit.2d
+ellipt->output->ParaView format: 1
+ellipt->output->ARH format: 1
+ellipt->output->AMDiS format: 1

extensions/demo/other/init/ellipt_krylov_precon.dat.2d

+dimension of world:             2
+
+elliptMesh->macro file name:    ./macro/macro.stand.2d
+elliptMesh->global refinements: 1	
+
+ellipt->mesh:                   elliptMesh
+ellipt->dim:                    2
+ellipt->components:             1
+ellipt->polynomial degree[0]:   1
+
+ellipt->solver:          fgmres
+ellipt->solver->max iteration:  100
+ellipt->solver->tolerance:      1.e-8
+ellipt->solver->info:           1
+ellipt->solver->print cycle:           1
+ellipt->solver->left precon:    solver
+ellipt->solver->left precon->solver: gmres
+ellipt->solver->left precon->solver->max iteration: 20
+ellipt->solver->left precon->solver->tolerance:      1.e-8
+ellipt->solver->left precon->solver->left precon: diag
+ellipt->solver->left precon->solver->restart: 30
+ellipt->solver->left precon->solver->info: 0
+ellipt->solver->right precon:   no
+
+ellipt->estimator[0]:              residual
+ellipt->estimator[0]->error norm:  1   % 1: H1_NORM, 2: L2_NORM
+ellipt->estimator[0]->C0:          0.1 % constant of element residual
+ellipt->estimator[0]->C1:          0.1 % constant of jump residual
+
+ellipt->marker[0]->strategy:       3   % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+
+ellipt->adapt[0]->tolerance:          1e-4
+ellipt->adapt->max iteration:         15
+
+ellipt->output->filename:       output/ellipt_mtl.2d
+%ellipt->output->filename:       output/ellipt_petsc.2d
+ellipt->output->ParaView format: 1

extensions/demo/other/init/parametric2.dat.2d

+dimension of world:      2
+
+parametricMesh->macro file name:       ./macro/macro.stand.2d
+parametricMesh->global refinements: 10    
+
+parametric->a: 0.2
+parametric->b: 1.0
+parametric->s: 1.3
+parametric->theta: M_PI/2
+
+parametric->mesh:                   parametricMesh
+parametric->dim:                    2
+parametric->polynomial degree[0]:   1
+parametric->components:             1
+
+parametric->solver:                cg
+parametric->solver->max iteration: 100
+parametric->solver->tolerance:     1.e-8
+parametric->solver->info:          2
+parametric->solver->left precon:   diag
+parametric->solver->right precon:  no
+
+parametric->estimator[0]:             residual % residual, recovery
+parametric->estimator[0]->error norm: 1   % 1: H1_NORM, 2: L2_NORM
+parametric->estimator[0]->C1:         0.1 % constant of jump residual
+
+parametric->marker[0]->strategy:      2   % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
+parametric->marker[0]->MSGamma:       0.5
+
+parametric->adapt[0]->tolerance:         1e-8
+parametric->adapt[0]->refine bisections: 1
+parametric->adapt[0]->info:              8
+
+parametric->adapt->max iteration:        8
+
+parametric->output->filename:         output/parametric2.2d
+parametric->output->ParaView format:  1
+parametric->output->append index:     0
+parametric->output->index length:     6
+parametric->output->index decimals:   3
+parametric->output->ARH format:  1
+parametric->output->AMDiS format: 1
+
+WAIT: 1
+

extensions/demo/other/src/cahnHilliard_navierStokes.cc

@@ -22,14 +22,14 @@
   #include "CouplingBaseProblem2.h"
 #endif
 #include "CahnHilliard.h"
-#include "ProblemStatMassConserve.h"
+#include "ProblemStatMassConserve2.h"
 #include "NavierStokes_TaylorHood.h"
 #include "Refinement.h"
 #include "MeshFunction_Level.h"
 
 using namespace AMDiS;
 
-typedef ::details::CahnHilliard<ProblemStatMassConserve> CHType;
+typedef ::detail::CahnHilliard<ProblemStatMassConserve2> CHType;
 typedef NavierStokes_TaylorHood NSType;
 
 class CahnHilliardNavierStokes : public CouplingBaseProblem< ProblemStat, CHType, NSType >
@@ -48,7 +48,8 @@ public: // methods
     chProb(&chProb_),
     nsProb(&nsProb_),
     refFunction(NULL),
-    refinement(NULL)
+    refinement(NULL),
+    minus1(-1.0)
   { }
 
   /// Destructor
@@ -67,16 +68,9 @@ public: // methods
 					chProb->getSolution()->getDOFVector(0));
   }
 
-  /// Called at the end of each timestep.
-  virtual void closeTimestep(AdaptInfo *adaptInfo)
-  {
-    super::closeTimestep(adaptInfo);
-    refinement->refine(1);
-  }
-
   /// Set initial condition and perform initial refinement
   virtual void solveInitialProblem(AdaptInfo *adaptInfo)
-  {
+  {    
     // initial refinement
     refinement->refine(0);
     
@@ -87,6 +81,35 @@ public: // methods
     }
 
     super::solveInitialProblem(adaptInfo);
+    
+    std::string prefix = "";
+    Parameters::get("prefix", prefix);
+    
+    double c = chProb->getSolution()->getDOFVector(0)->Int();
+    std::ofstream out(("concentration" + prefix + ".dat").c_str(), ios_base::out);
+    out.setf(std::ios::scientific, std::ios::floatfield );
+    out.precision(16);
+    out << c << "\n";
+    
+    // activate mass-conservative coarsening
+    refinement->setRefineOperation(adaptInfo, chProb->getProblem());
+  }
+
+  /// Called at the end of each timestep.
+  virtual void closeTimestep(AdaptInfo *adaptInfo)
+  {
+    super::closeTimestep(adaptInfo);
+    
+    std::string prefix = "";
+    Parameters::get("prefix", prefix);
+    
+    double c = chProb->getSolution()->getDOFVector(0)->Int();
+    std::ofstream out(("concentration" + prefix + ".dat").c_str(), ios_base::app);
+    out.setf(std::ios::scientific, std::ios::floatfield );
+    out.precision(16);
+    out << c << "\n";
+    
+    refinement->refine(1);
   }
   
   virtual void fillCouplingOperators()
@@ -105,13 +128,14 @@ public: // methods
       nsProb->getProblem(0)->addVectorOperator(opCGradMu, i);
     }
     
-    // < u * grad(c) , theta >
+    // < - u * c , grad(theta) >
     for (size_t i = 0; i < dow; i++) {
       Operator *opUGradC = new Operator(chProb->getFeSpace(1), 
-					chProb->getFeSpace(0));
+					chProb->getFeSpace(0));      
       opUGradC->addTerm(new VecAndPartialDerivative_FOT(
-	nsProb->getSolution()->getDOFVector(i), i), GRD_PHI);
-      chProb->getProblem()->addMatrixOperator(opUGradC, 1, 0);
+	nsProb->getSolution()->getDOFVector(i), // u_i
+	i), GRD_PSI);
+      chProb->getProblem()->addMatrixOperator(opUGradC, 1, 0, &minus1, &minus1);
     }
   }
 
@@ -121,6 +145,9 @@ protected: // variables
   NSType *nsProb;
   PhaseFieldRefinement* refFunction;
   RefinementLevelDOF *refinement;
+
+  double minus1;
+  double concentration;
 };