some new baseproblems

extensions/base_problems/BaseProblem.h

@@ -14,7 +14,7 @@ using namespace AMDiS;
 const Flag MESH_ADOPTED = 1<<2;
 const Flag DATA_ADOPTED = 1<<3;
 
-template<typename ProblemType=ProblemStat, bool safeguard = false>
+template<class ProblemType=ProblemStat, bool safeguard = false>
 class BaseProblem : public ProblemIterationInterface,
 		    public ProblemInstatBase
 {
@@ -124,12 +124,19 @@ public:
     return NULL;
   }
 
-  ProblemType *getProblem(std::string name)
+  ProblemType *getProblem(std::string name_)
   {
-    if (name == "prob")
+    if (name_ == "prob")
       return prob;
     else
-      throw(std::runtime_error("problem with given name '" + name + "' does not exist"));
+      throw(std::runtime_error("problem with given name '" + name_ + "' does not exist"));
+  }
+
+  void setProblem(ProblemType *prob_, bool destroy = true)
+  {
+    if (destroy && prob)
+      delete prob;
+    prob = prob_;
   }
 
   // setting methods
@@ -168,6 +175,7 @@ public:
   static void initFileWriterFromFile(AdaptInfo* adaptInfo, FileWriterTemplated<T>& fileWriter, bool keep_all = false);
   
 protected:
+  std::string name;
 
   ProblemType *prob;
 
@@ -183,7 +191,6 @@ protected:
   double oldTimestep;
 
 private:
-
   std::vector<std::pair<int,int> > fixedMatrixTimestep;
 
 };

extensions/base_problems/BaseProblem.hh

@@ -3,6 +3,7 @@ using namespace AMDiS;
 template<typename ProblemType, bool safeguard>
 BaseProblem<ProblemType, safeguard>::BaseProblem(const std::string &name_, bool createProblem) 
   : ProblemInstatBase(name_,NULL),
+    name(name_),
     prob(NULL),
     dim(1),
     dow(1),
@@ -12,7 +13,7 @@ BaseProblem<ProblemType, safeguard>::BaseProblem(const std::string &name_, bool
 {
   // create basic problems
   if (createProblem)
-    prob = new ProblemType(name + "->space");
+    prob = new ProblemType(name_ + "->space");
   
   dow = Global::getGeo(WORLD);
 }
@@ -27,6 +28,7 @@ void BaseProblem<ProblemType, safeguard>::initialize(Flag initFlag,
   TEST_EXIT(prob)("Problem must be created! Add flag createProblem=true to constructor!\n");
   
   prob->initialize(initFlag, adoptProblem, adoptFlag);
+  ProblemInstatBase::initialize(INIT_ALL);
   dim = getMesh()->getDim();
 }
 

extensions/base_problems/BaseProblem_RB.h

@@ -23,19 +23,57 @@ public: // typedefs
 
 public:
   
-  BaseProblem_RB(std::string name_) : super(name_), minus1(-1.0) {}
+  BaseProblem_RB(std::string name_, bool createProblem = true) : super(name_, createProblem), minus1(-1.0) {}
   
-  void acceptTimestep() { prob->acceptTimestep(); }
-  void reset() { prob->reset(); }
+  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) 
+  {
+    prob->oneIteration(adaptInfo, toDo);
+    
+    oldTimestep = adaptInfo->getTimestep();
+    oldMeshChangeIdx = getMesh()->getChangeIndex();
+  }
   
-  void setRosenbrockMethod(RosenbrockMethod *method) { prob->setRosenbrockMethod(method); }
+  Flag stageIteration(AdaptInfo *adaptInfo, Flag flag, bool asmMatrix, bool asmVector)
+  {
+    return prob->stageIteration(adaptInfo, flag, asmMatrix, asmVector);
+  }
+    
+  double estimateTimeError(AdaptInfo* adaptInfo)
+  {
+    return prob->estimateTimeError(adaptInfo);
+  }
+  
+  void acceptTimestep(AdaptInfo* adaptInfo) { prob->acceptTimestep(adaptInfo); }  
+    
+  double getNewTimestep(double tol, double tau, bool restrict = true)
+  {
+    return prob->getNewTimestep(tol, tau, restrict);
+  }
+    
+  double getNewTimestep(double tol, double oldTau, double tau, bool restrict = true)
+  {
+    return prob->getNewTimestep(tol, oldTau, tau, restrict);
+  }
+  
+  void reduceOrder(double tauRej, double tau) { prob->reduceOrder(tauRej, tau); }
   void setTau(double *ptr) { prob->setTau(ptr); }
-  void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3, double *ptr4) { prob->setTauGamma(ptr0,ptr1,ptr2,ptr3,ptr4); }
-
-  double* getTauGamma() { return prob->getTauGamma(); }
-  double* getMinusTauGamma() { return prob->getMinusTauGamma(); }
-  double* getInvTauGamma() { return prob->getInvTauGamma(); }
-  double* getMinusInvTauGamma() { return prob->getMinusInvTauGamma(); }
+  void setOldTime(double t) { prob->setOldTime(t); }
+    
+  /// get the i-th stage-solution-vector
+  DOFVector<double> *getStageSolution(int i)
+  { FUNCNAME("CoupledBaseProblem::getStageSolution(i,j)");
+  
+    TEST_EXIT(0<=i && i<getNumComponents())("Indices out of range!\n");
+    return prob->getStageSolution(i);
+  }  
+  
+  /// get the i-th unVec-vector
+  DOFVector<double> *getUnVec(int i) 
+  { FUNCNAME("CoupledBaseProblem::getUnVec(i)");
+  
+    TEST_EXIT(0<=i && i<getNumComponents())("Indices out of range!\n");
+    return prob->getUnVec(i);
+  }
   
 protected:
   double minus1;

extensions/base_problems/BaseProblem_RB0.h

+/** \file BaseProblem_RB0.h */
+
+#ifndef BASE_PROBLEM_RB0_H
+#define BASE_PROBLEM_RB0_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+// #include "time/ExtendedRosenbrockStationary.h"
+
+using namespace AMDiS;
+
+/** BaseProblems for Rosenbrock time-discretization
+ */
+
+class BaseProblem_RB0 : public BaseProblem<RosenbrockStationary>
+{
+public: // typedefs
+
+  typedef BaseProblem<RosenbrockStationary> super;
+
+// Rosenbrock methods
+//_________________________________________________________________________
+
+public:
+  
+  BaseProblem_RB0(std::string name_, bool createProblem = true) : super(name_, createProblem), minus1(-1.0) {}
+  
+  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo) 
+  { 
+    prob->oneIteration(adaptInfo, toDo);
+    
+    oldTimestep = adaptInfo->getTimestep();
+    oldMeshChangeIdx = getMesh()->getChangeIndex();
+  }
+  
+  void setRosenbrockMethod(RosenbrockMethod *method) { prob->setRosenbrockMethod(method); }
+  void setTau(double *ptr) { prob->setTau(ptr); }
+  void setTauGamma(double *ptr0, double *ptr1, double *ptr2, double *ptr3)
+  {
+    prob->setTauGamma(ptr0, ptr1, ptr2, ptr3);
+  }  
+    
+  /// get the i-th stage-solution-vector
+  DOFVector<double> *getStageSolution(int i)
+  { FUNCNAME("CoupledBaseProblem::getStageSolution(i,j)");
+  
+    TEST_EXIT(0<=i && i<getNumComponents())("Indices out of range!\n");
+    return prob->getStageSolution(i);
+  }  
+  
+  /// get the i-th unVec-vector
+  DOFVector<double> *getUnVec(int i) 
+  { FUNCNAME("CoupledBaseProblem::getUnVec(i)");
+  
+    TEST_EXIT(0<=i && i<getNumComponents())("Indices out of range!\n");
+    return prob->getUnVec(i);
+  }
+  
+protected:
+  double minus1;
+};
+
+#endif // BASE_PROBLEM_RB_H
\ No newline at end of file

extensions/base_problems/CahnHilliardNavierStokes_RB.cc

@@ -272,24 +272,24 @@ void CahnHilliardNavierStokes_RB::fillOperators()
   // -----------------------------------------------------------------------------------
   for (size_t i = 0; i < dow; ++i) {
     /// < d_t(u_i) , psi >
-    prob->addTimeOperator(2+i, 2+i);
+    prob->addTimeOperator(2+i, 2+i, density);
  
     /// < u_s*grad(u_s_i) , psi >
     Operator *opUGradU = new Operator(getFeSpace(2+i), getFeSpace(2+i));
-    opUGradU->addTerm(new WorldVec_FOT(stage_velocity, -1.0), GRD_PHI);
+    opUGradU->addTerm(new WorldVec_FOT(stage_velocity, -density), GRD_PHI);
     opUGradU->setUhOld(stage_velocity[i]);
     prob->addVectorOperator(*opUGradU, 2+i);
 
     /// < u_old*grad(u_i) , psi >
     Operator *opUGradV = new Operator(getFeSpace(2+i), getFeSpace(2+i));
-    opUGradV->addTerm(new WorldVec_FOT(un_velocity, -1.0), GRD_PHI);
-    prob->addMatrixOperator(*opUGradV, 2+i, 2+i, &minus1, &minus1);
+    opUGradV->addTerm(new WorldVec_FOT(un_velocity, density), GRD_PHI);
+    prob->addMatrixOperator(*opUGradV, 2+i, 2+i);
 
     /// < u*grad(u_old_i) , psi >
     for (size_t j = 0; j < dow; ++j) {
       Operator *opVGradU = new Operator(getFeSpace(2+i), getFeSpace(2+j));
-      opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, -1.0));
-      prob->addMatrixOperator(*opVGradU, 2+i, 2+j, &minus1, &minus1);
+      opVGradU->addTerm(new PartialDerivative_ZOT(un_velocity[i], j, density));
+      prob->addMatrixOperator(*opVGradU, 2+i, 2+j);
     }
 
     /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
@@ -297,15 +297,15 @@ void CahnHilliardNavierStokes_RB::fillOperators()
   
     /// < p , d_i(psi) >
     Operator *opGradP = new Operator(getFeSpace(2+i),getFeSpace(2+dow));
-    opGradP->addTerm(new PartialDerivative_FOT(i), GRD_PSI);
+    opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
     opGradP->setUhOld(prob->getStageSolution(2+dow));
     prob->addVectorOperator(*opGradP, 2+i);
     prob->addMatrixOperator(*opGradP, 2+i, 2+dow, &minus1, &minus1);
   
     /// external force, i.e. gravitational force
-    if (norm(force) > DBL_TOL) {
+    if (std::abs(force[i]) > DBL_TOL) {
       Operator *opForce = new Operator(getFeSpace(2+i), getFeSpace(2+i));
-      opForce->addTerm(new Simple_ZOT(force[i]));
+      opForce->addTerm(new Simple_ZOT(density*force[i]));
       prob->addVectorOperator(*opForce, 2+i);
     }
     
@@ -343,6 +343,10 @@ void CahnHilliardNavierStokes_RB::fillOperators()
     opDivU->setUhOld(stage_velocity[i]);
     prob->addVectorOperator(*opDivU, 2+dow);
   }
+  
+  Operator *opZero = new Operator(getFeSpace(2+dow),getFeSpace(2+dow));
+  opZero->addTerm(new Simple_ZOT(0.0));
+  prob->addMatrixOperator(*opZero, 2+dow, 2+dow);
 }
 
 
@@ -353,7 +357,7 @@ void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
     if (laplaceType == 1) {
       for (unsigned j = 0; j < dow; ++j) {
         Operator *opLaplaceUi1 = new Operator(getFeSpace(2+i), getFeSpace(2+j));
-        opLaplaceUi1->addTerm(new MatrixIJ_SOT(1-i, 1-j, viscosity));
+        opLaplaceUi1->addTerm(new MatrixIJ_SOT(j, i, viscosity));
         prob->addMatrixOperator(*opLaplaceUi1, 2+i, 2+j);
 	opLaplaceUi1->setUhOld(prob->getStageSolution(2+j));
 	prob->addVectorOperator(*opLaplaceUi1, 2+i, &minus1, &minus1);
@@ -362,7 +366,7 @@ void CahnHilliardNavierStokes_RB::addLaplaceTerm(int i)
     
     /// < nu*grad(u'_i) , grad(psi) >
     Operator *opLaplaceUi = new Operator(getFeSpace(2+i), getFeSpace(2+i));
-      opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+    opLaplaceUi->addTerm(new Simple_SOT(viscosity));
     prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
     opLaplaceUi->setUhOld(prob->getStageSolution(2+i));
     prob->addVectorOperator(*opLaplaceUi, 2+i, &minus1, &minus1);

extensions/base_problems/CoupledBaseProblem.h

+/** \file CoupledBaseProblem.h */
+
+#ifndef COUPLED_BASE_PROBLEM_H
+#define COUPLED_BASE_PROBLEM_H
+
+#include "AMDiS.h"
+
+// coupling structures
+#include "CouplingIterationInterface.h"
+#include "CouplingTimeInterface.h"
+#include "CouplingProblemStat.h"
+#include "BaseProblem.h"
+#include "ProblemInstat.h"
+
+using namespace AMDiS;
+
+
+/**
+  * \ingroup Problem 
+  *
+  * \brief
+  */
+template<class ProblemType=ProblemStat, class BaseProblemType=BaseProblem<ProblemStat> >
+class CoupledBaseProblem : public CouplingIterationInterface,
+			   public CouplingTimeInterface,
+			   public CouplingProblemStatImpl<ProblemType>
+{
+public:
+  typedef CouplingProblemStatImpl<ProblemType> CProblemStat;
+
+  CoupledBaseProblem(std::string name_, 
+		     BaseProblemType *prob0_,
+		     BaseProblemType *prob1_=NULL, 
+		     BaseProblemType *prob2_=NULL, 
+		     BaseProblemType *prob3_=NULL, 
+		     BaseProblemType *prob4_=NULL)
+  : CProblemStat(name_)
+  {
+    dow = Global::getGeo(WORLD);
+    baseProblems.push_back(prob0_);
+    if (prob1_) baseProblems.push_back(prob1_);
+    if (prob2_) baseProblems.push_back(prob2_);
+    if (prob3_) baseProblems.push_back(prob3_);
+    if (prob4_) baseProblems.push_back(prob4_);
+  }
+
+  ~CoupledBaseProblem() { }  
+  
+  virtual void initialize(Flag initFlag,
+			  ProblemStat *adoptProblem = NULL,
+			  Flag adoptFlag = INIT_NOTHING)
+  {  
+    for (size_t i = 0; i < baseProblems.size(); i++) {
+      for (size_t j = 0; j < baseProblems[i]->getNumProblems(); j++)
+	addProblem(dynamic_cast<ProblemType*>(baseProblems[i]->getProblem(j)));
+      
+      addIterationInterface(baseProblems[i]);
+      addTimeInterface(baseProblems[i]);
+    }
+
+    // initialize all ProblemStat
+    CProblemStat::initialize(initFlag, adoptProblem, adoptFlag);
+    
+    // set mesh dimension
+    dim = CProblemStat::getMesh(0)->getDim();
+  }  
+  
+  virtual void initData() {}  
+  
+  virtual void initTimeInterface()
+  {
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->initData();
+    initData();
+    
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->fillOperators();
+    fillCouplingOperators();
+    
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->fillBoundaryConditions();
+    fillCouplingBoundaryConditions();
+  }
+  
+  
+  virtual void fillCouplingOperators() {}
+  virtual void fillCouplingBoundaryConditions() {}
+  
+  
+  /// get the j-th solution-vector of the i-th problem
+  DOFVector<double> *getSolution(int i, int j)
+  { FUNCNAME("CoupledBaseProblem::getSolution(i,j)");
+    TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
+      ("Indices out of range!\n");
+    return baseProblems[i]->getSolution()->getDOFVector(j);
+  }
+  
+  
+  /// pointer to the j-th feSpace of the i-th problem
+  inline const FiniteElemSpace* getFeSpace(int i, int j=0) 
+  { FUNCNAME("CoupledBaseProblem::getFeSpace(i,j)");
+    TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
+      ("Indices out of range!\n");
+    return baseProblems[i]->getFeSpace(j); 
+  }
+
+protected:
+  std::vector<BaseProblemType*> baseProblems;
+
+  unsigned dim;	// dimension of the meshes
+  unsigned dow;			// dimension of the world
+};
+
+
+#endif // COUPLED_BASE_PROBLEM_H

extensions/base_problems/CoupledBaseProblem_RB.h

+/** \file CoupledBaseProblem_RB.h */
+
+#ifndef COUPLED_BASE_PROBLEM_RB_H
+#define COUPLED_BASE_PROBLEM_RB_H
+
+#include "AMDiS.h"
+
+// coupling structures
+#include "CouplingIterationInterface.h"
+#include "CouplingTimeInterface.h"
+#include "CouplingProblemStat.h"
+#include "BaseProblem_RB.h"
+#include "CoupledBaseProblem.h"
+#include "time/ExtendedRosenbrockStationary.h"
+
+using namespace AMDiS;
+
+
+/**
+  * \ingroup Problem 
+  *
+  * \brief
+  */
+class CoupledBaseProblem_RB : public CoupledBaseProblem<ExtendedRosenbrockStationary, BaseProblem_RB>
+{
+public:
+  typedef CoupledBaseProblem<ExtendedRosenbrockStationary, BaseProblem_RB> super;
+  typedef BaseProblem_RB BaseProblemType;
+  typedef CouplingProblemStatImpl<ExtendedRosenbrockStationary> CProblemStat;
+
+  CoupledBaseProblem_RB(std::string name_, 
+			BaseProblemType *prob0_,
+			BaseProblemType *prob1_=NULL, 
+			BaseProblemType *prob2_=NULL, 
+			BaseProblemType *prob3_=NULL, 
+			BaseProblemType *prob4_=NULL)
+  : super(name_,prob0_,prob1_,prob2_,prob3_,prob4_)
+  {
+    int componentShift = 0;
+    for (size_t i = 0; i < baseProblems.size(); i++) {
+      baseProblems[i]->setProblem(new ExtendedRosenbrockStationary(baseProblems[i]->getName() + "->space", componentShift));
+      componentShift += baseProblems[i]->getNumComponents();
+    }
+  }
+
+  ~CoupledBaseProblem_RB() { }  
+  
+  /// get the j-th stage-solution-vector of the i-th problem
+  DOFVector<double> *getStageSolution(int i, int j)
+  { FUNCNAME("CoupledBaseProblem::getStageSolution(i,j)");
+  
+    TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
+      ("Indices out of range!\n");
+    return baseProblems[i]->getStageSolution(j);
+  }  
+  
+  /// get the j-th unVec-vector of the i-th problem
+  DOFVector<double> *getUnVec(int i, int j=0) 
+  { FUNCNAME("CoupledBaseProblem::getUnVec(i,j)");
+  
+    TEST_EXIT(0<=i && 0<=j && i<baseProblems.size() && j<=baseProblems[i]->getNumComponents())
+      ("Indices out of range!\n");
+    return baseProblems[i]->getUnVec(j);
+  }
+  
+  virtual double estimateTimeError(AdaptInfo* adaptInfo)
+  {
+    double errorEst = 0.0;
+    for (size_t i = 0; i < CProblemStat::getNumProblems(); i++)
+      errorEst = std::max(errorEst, CProblemStat::getProblem(i)->estimateTimeError(adaptInfo));
+    
+    adaptInfo->setTimeEst(errorEst);
+    return errorEst;
+  }
+  
+  double getNewTimestep(double tol, double tau, bool restrict = true)
+  {
+    double newTimestep = 1.e15;
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      newTimestep = std::min(newTimestep, baseProblems[i]->getNewTimestep(tol, tau, restrict));
+    return newTimestep;
+  }
+  
+  double getNewTimestep(double tol, double oldTau, double tau, bool restrict = true)
+  {
+    double newTimestep = 1.e15;
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      newTimestep = std::min(newTimestep, baseProblems[i]->getNewTimestep(tol, oldTau, tau, restrict));
+    return newTimestep;
+  }
+
+  // Rosenbrock methods
+  //_________________________________________________________________________
+
+  void acceptTimestep(AdaptInfo* adaptInfo)
+  {
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->acceptTimestep(adaptInfo);
+  }
+
+//   void reset()
+//   {
+//     for (size_t i = 0; i < CProblemStat::getNumProblems(); i++)
+//       CProblemStat::getProblem(i)->reset();
+//     for (size_t i = 0; i < baseProblems.size(); i++)
+//       baseProblems[i]->reset();
+//   }
+
+  void setTau(double *ptr)
+  {
+    for (size_t i = 0; i < CProblemStat::getNumProblems(); i++)
+      CProblemStat::getProblem(i)->setTau(ptr);
+  }
+  
+  void setOldTime(double t)
+  {
+    for (size_t i = 0; i < CProblemStat::getNumProblems(); i++)
+      CProblemStat::getProblem(i)->setOldTime(t);
+  }
+  
+  void setTime(AdaptInfo* adaptInfo)
+  {
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->setTime(adaptInfo);
+  }
+	
+  void reduceOrder(double tauRej, double tau)
+  {
+    for (size_t i = 0; i < baseProblems.size(); i++)
+      baseProblems[i]->reduceOrder(tauRej, tau);
+  }
+  
+  double* getOldTimePtr()
+  {
+    return CProblemStat::getProblem(0)->getOldTimePtr();
+  }
+};
+
+
+#endif // COUPLED_BASE_PROBLEM_RB_H

extensions/base_problems/NavierStokes_TH_MultiPhase.cc

 #include "NavierStokes_TH_MultiPhase.h"
 
-using namespace std;
 using namespace AMDiS;
 
-NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_) :
-  super(name_),
+NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_, bool createProblem) :
+  super(name_, createProblem),
   multiPhase(NULL),
   densityPhase(NULL),
   viscosityPhase(NULL),
@@ -12,7 +11,8 @@ NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_)
   viscosity2(1.0),
   density1(1.0),
   density2(1.0)
-{
+{ FUNCNAME("NavierStokes_TH_MultiPhase::_constructor()");
+
   Initfile::get(name + "->viscosity1", viscosity1); // viscosity of fluid 1
   Initfile::get(name + "->viscosity2", viscosity2); // viscosity of fluid 2
   Initfile::get(name + "->density1", density1); // density of fluid 1
@@ -22,35 +22,31 @@ NavierStokes_TH_MultiPhase::NavierStokes_TH_MultiPhase(const std::string &name_)
     viscosity1 = viscosity;
   
   TEST_EXIT(viscosity1 > 0.0)("Positive viscosity is necessary!\n");
-};
+}
 
 
 NavierStokes_TH_MultiPhase::~NavierStokes_TH_MultiPhase()
-{ FUNCNAME("NavierStokes_TH_MultiPhase::~NavierStokes_TH_MultiPhase()");
-
-  //delete multiPhase;
-  delete densityPhase;
-  delete viscosityPhase;
-};
+{
+  if (multiPhase) { delete multiPhase; multiPhase = NULL; }
+  if (densityPhase) { delete densityPhase; densityPhase = NULL; }
+  if (viscosityPhase) { delete viscosityPhase; viscosityPhase = NULL; }
+}
 
 
 void NavierStokes_TH_MultiPhase::initData()