ProblemStatMassConserve.h 8.55 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
/******************************************************************************
 *
 * 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 ProblemStatMassConserve.h */


#ifndef PROBLEM_STAT_MASS_CONSERVE
#define PROBLEM_STAT_MASS_CONSERVE

#include "ExtendedProblemStat.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 ProblemStatMassConserve : public ExtendedProblemStat
{
  /// constructor: create a temporary problem that is initialized with 
  /// similar parameters as the regular one, except one additional component
  /// with an extra mesh for a temporary oldSolution
  ProblemStatMassConserve(std::string name_)
  : ExtendedProblemStat(name_),
    initialMeshAdoption(false),
    prob2(nullptr),
    comp(0)
  {
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    TEST_EXIT(false)("Multi-Mesh approach does not work in parallel mode for the moment. Need to be implemented!\n");
#endif
    
    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", two);
    Parameters::set(name2 + "->polynomial degree[0]", degree);
    Parameters::set(name2 + "->polynomial degree[1]", degree);
    Parameters::set(name2 + "->dim", dim);
    Parameters::set(name2 + "->mesh", meshName);
    Parameters::set(name2 + "->refinement set[0]", one);
    Parameters::set(name2 + "->refinement set[1]", two);
#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);
  }
  
  
  ~ProblemStatMassConserve()
  {
    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));
      
    int refSet = -1;
    Parameters::get(name + "_tmp->refinement set[0]", refSet);
    if (refSet < 0)
      refSet = 0;
    
    std::map<int, Mesh*> meshForRefinementSet;
    meshForRefinementSet[refSet] = meshes2[0];
        
    for (int i = 1; i < prob2->getNumComponents(); i++) {
      refSet = -1;
      Parameters::get(name + "_tmp->refinement set[" + lexical_cast<string>(i) + "]", 
		      refSet);
      if (refSet < 0)
	refSet = 0;
      
      if (meshForRefinementSet[refSet] == nullptr) {
	Mesh *newMesh = new Mesh(getMesh(comp)->getName(), getMesh(comp)->getDim());
	meshForRefinementSet[refSet] = newMesh;
	meshes2.push_back(newMesh);
      }
      prob2->setComponentMesh(i, meshForRefinementSet[refSet]);
    }
    prob2->setMeshes(meshes2);    
    prob2->initialize(INIT_ALL - INIT_MESH);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
128
    for (size_t i = 0; i < meshes2.size(); i++) {
129 130 131 132 133 134 135 136 137 138 139 140 141
      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.
Praetorius, Simon's avatar
Praetorius, Simon committed
142
      if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) {
143 144
	meshes2[i]->initialize();
	prob2->getRefinementManager()->globalRefine(meshes2[i], globalRefinements);
Praetorius, Simon's avatar
Praetorius, Simon committed
145
      }
146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
    }
    
    
    for (int i = 0; i < prob2->getNumComponents(); i++)
      prob2->getSolution(i)->setCoarsenOperation(NO_OPERATION);
    
    fillOperators();
  }
  
  
  /// in the first iteration copy initial-solution to temporary problem
  /// and prepare the temporary mesh
  void buildBeforeRefine(AdaptInfo* adaptInfo, 
			 Flag markFlag) override
  {
    if (!initialMeshAdoption) {
      MeshStructure meshStructure2;
      meshStructure2.init(getMesh(comp));
      prob2->getCoarseningManager()->globalCoarsen(prob2->getMesh(1), -2);
      meshStructure2.fitMeshToStructure(prob2->getMesh(1), prob2->getRefinementManager());
      
      prob2->getSolution(1)->interpol(getSolution(comp));
      initialMeshAdoption = true;
    }
    
    ExtendedProblemStat::buildBeforeRefine(adaptInfo, markFlag);
  }
  
  
  /// 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
  {  
    prob2->buildAfterCoarsen(adaptInfo, markFlag, true, true);
    solve2(adaptInfo);
    getSolution(comp)->copy(*prob2->getSolution(0));
    
    // call original buildAfterCoarsen
    ExtendedProblemStat::buildAfterCoarsen(adaptInfo, markFlag, asmM, asmV); 
  }
  
  
  /// copy solution to temporary problem after regular solve
  void solve(AdaptInfo *adaptInfo, 
	     bool createMatrixData = true, 
	     bool storeMatrixData = false) override
  {
    ExtendedProblemStat::solve(adaptInfo, createMatrixData, storeMatrixData);
    
    MeshStructure meshStructure2;
    meshStructure2.init(getMesh(comp));
    prob2->getCoarseningManager()->globalCoarsen(prob2->getMesh(1), -2);
    meshStructure2.fitMeshToStructure(prob2->getMesh(1), prob2->getRefinementManager());
    
    prob2->getSolution(1)->interpol(getSolution(comp));
  }
  
  /// solve only one block of temporary problem
  void solve2(AdaptInfo *adaptInfo, 
	      bool createMatrixData = true, 
	      bool storeMatrixData = false)
  {
    SolverMatrix<Matrix<DOFMatrix*> > solverMatrix_;
    Matrix<DOFMatrix*> mat(1,1);
    mat[0][0] = prob2->getSystemMatrix(0,0);
    solverMatrix_.setMatrix(mat);
    
    std::vector<const FiniteElemSpace*> feSpaces0; 
    feSpaces0.push_back(prob2->getFeSpace(0));
    
    SystemVector vec_sol("solution_0", feSpaces, 1, false);
    vec_sol.setDOFVector(0, prob2->getSolution(0));
    SystemVector vec_rhs("rhs_0", feSpaces, 1, false);
    vec_rhs.setDOFVector(0, prob2->getRhsVector(0));
    
    // solve only one block of the system
    prob2->getSolver()->solveSystem(solverMatrix_, vec_sol, vec_rhs, 
			createMatrixData, storeMatrixData);
  }

  
  /// add mass-matrix and transfer-operator to temporary problem, to
  /// project solution to new mesh
  void fillOperators()
  {
    // conversion propblem 2
    const FiniteElemSpace* feSpace2_0 = prob2->getFeSpace(0);
    const FiniteElemSpace* feSpace2_1 = prob2->getFeSpace(1);

    Operator *opMnew0 = new Operator(feSpace2_0, feSpace2_0);
    opMnew0->addTerm(new Simple_ZOT);
    Operator *opMold01 = new Operator(feSpace2_0, feSpace2_1);
    opMold01->addTerm(new VecAtQP_ZOT(prob2->getSolution(1)));
    
    // assemble only one block of the system
    prob2->addMatrixOperator(*opMnew0,0,0);
    prob2->addVectorOperator(*opMold01, 0); 
  }
  
private:
  /// true if initial-solution/-mesh is copied to temporary problem
  bool initialMeshAdoption; 	
  
  /// 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