ProblemStatMassConserve.h 8.55 KB
 Praetorius, Simon committed Oct 23, 2013 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(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 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 meshForRefinementSet; meshForRefinementSet[refSet] = meshes2[0]; for (int i = 1; i < prob2->getNumComponents(); i++) { refSet = -1; Parameters::get(name + "_tmp->refinement set[" + lexical_cast(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 committed Oct 29, 2013 128  for (size_t i = 0; i < meshes2.size(); i++) {  Praetorius, Simon committed Oct 23, 2013 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 committed Oct 29, 2013 142  if (initMesh && meshes2[i] && !(meshes2[i]->isInitialized())) {  Praetorius, Simon committed Oct 23, 2013 143 144  meshes2[i]->initialize(); prob2->getRefinementManager()->globalRefine(meshes2[i], globalRefinements);  Praetorius, Simon committed Oct 29, 2013 145  }  Praetorius, Simon committed Oct 23, 2013 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 > solverMatrix_; Matrix mat(1,1); mat[0][0] = prob2->getSystemMatrix(0,0); solverMatrix_.setMatrix(mat); std::vector 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