Base Problems added

extensions/base_problems/BaseProblem_RB.h

+/** \file BaseProblem_RB.h */
+
+#ifndef BASE_PROBLEM_RB_H
+#define BASE_PROBLEM_RB_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "time/ExtendedRosenbrockStationary.h"
+
+using namespace AMDiS;
+
+/** BaseProblems for Rosenbrock time-discretization
+ */
+
+class BaseProblem_RB : public BaseProblem<ExtendedRosenbrockStationary>
+{
+public: // typedefs
+
+  typedef BaseProblem<ExtendedRosenbrockStationary> super;
+
+// Rosenbrock methods
+//_________________________________________________________________________
+
+public:
+  
+  BaseProblem_RB(std::string name_) : super(name_), minus1(-1.0) {}
+  
+  void acceptTimestep() { prob->acceptTimestep(); }
+  void reset() { prob->reset(); }
+  
+  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); }
+
+  double* getTauGamma() { return prob->getTauGamma(); }
+  double* getMinusTauGamma() { return prob->getMinusTauGamma(); }
+  double* getInvTauGamma() { return prob->getInvTauGamma(); }
+  double* getMinusInvTauGamma() { return prob->getMinusInvTauGamma(); }
+  
+protected:
+  double minus1;
+};
+
+#endif // BASE_PROBLEM_RB_H
\ No newline at end of file

extensions/base_problems/CahnHilliardNavierStokes.cc

+#include "CahnHilliardNavierStokes.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "POperators.h"
+
+using namespace AMDiS;
+
+CahnHilliardNavierStokes::CahnHilliardNavierStokes(const std::string &name_) :
+  super(name_),
+  useMobility(false),
+  useConservationForm(false),
+  doubleWell(0),
+  gamma(1.0),
+  eps(0.1),
+  minusEps(-0.1),
+  epsInv(10.0),
+  minusEpsInv(-10.0),
+  epsSqr(0.01),
+  minusEpsSqr(-0.01),
+  laplaceType(0),
+  nonLinTerm(2),
+  oldTimestep(0.0),
+  viscosity(1.0),
+  density(1.0),
+  c(0.0),
+  sigma(1.0),
+  surfaceTension(1.0),
+  reinit(NULL),
+  velocity(NULL),
+  fileWriter(NULL)
+{
+  // parameters for CH
+  Parameters::get(name + "->use mobility", useMobility); // mobility
+  Parameters::get(name + "->gamma", gamma); // mobility
+  Parameters::get(name + "->epsilon", eps); // interface width
+  
+  // type of double well: 0= [0,1], 1= [-1,1]
+  Parameters::get(name + "->double-well type", doubleWell); 
+
+  // parameters for navier-stokes
+  Initfile::get(name + "->viscosity", viscosity);
+  Initfile::get(name + "->density", density);
+  // type of laplace operator: 0... div(nu*grad(u)), 1... div(0.5*nu*(grad(u)+grad(u)^T))
+  Initfile::get(name + "->laplace operator", laplaceType); 
+  // type of non-linear term: 0... u^old*grad(u_i^old), 1... u'*grad(u_i^old), 2... u^old*grad(u'_i)
+  Initfile::get(name + "->non-linear term", nonLinTerm); 
+
+  // Parameters for CH-Coupling
+  Initfile::get(name + "->use conservation form", useConservationForm);
+
+  // Parameters for NS-Coupling
+  Initfile::get(name + "->sigma", sigma); // surface tension
+  surfaceTension = sigma*3.0/(2.0*sqrt(2.0)) / eps;
+
+    
+  force.set(0.0);
+  Initfile::get(name + "->force", force);
+  
+  // transformation of the parameters
+  minusEps = -eps;
+  epsInv = 1.0/eps;
+  minusEpsInv = -epsInv;
+  epsSqr = sqr(eps);
+  minusEpsSqr = -epsSqr;
+}
+
+
+CahnHilliardNavierStokes::~CahnHilliardNavierStokes() 
+{ FUNCNAME("CahnHilliardNavierStokes::~CahnHilliardNavierStokes()");
+
+  if (reinit != NULL) {
+    delete reinit;
+    reinit = NULL;
+  }
+  if (velocity != NULL) {
+    delete velocity;
+    velocity = NULL;
+  }
+  delete fileWriter;
+  fileWriter = NULL;
+}
+
+
+void CahnHilliardNavierStokes::initData()
+{ FUNCNAME("CahnHilliardNavierStokes::initData()");
+
+  dim = getMesh()->getDim();
+  // create instance redistancing class
+  reinit = new HL_SignedDistTraverse("reinit", dim);
+  
+  if (velocity == NULL)
+    velocity = new DOFVector<WorldVector<double> >(getFeSpace(2), "velocity");
+  fileWriter = new FileVectorWriter(name + "->velocity->output", getFeSpace(2)->getMesh(), velocity);
+
+  super::initData();
+}
+
+
+void CahnHilliardNavierStokes::fillOperators()
+{ FUNCNAME("CahnHilliardNavierStokes::fillOperators()");
+  MSG("CahnHilliardNavierStokes::fillOperators()\n");
+  
+  // variable order:
+  // (c, mu, u0, u1 [, u2], p)
+  
+  WorldVector<DOFVector<double>* > vel;
+  for (size_t i = 0; i < dow; i++)
+    vel[i] = prob->getSolution()->getDOFVector(2+i);
+
+
+  // dt(c) = laplace(mu) - u*grad(c)
+  // -----------------------------------
+  /// < phi*c/tau , psi >
+  Operator *opChMnew = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+  opChMnew->addZeroOrderTerm(new Simple_ZOT);
+  prob->addMatrixOperator(*opChMnew, 0, 0, getInvTau()); 
+  
+  /// < phi*c^old/tau , psi >
+  Operator *opChMold = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+  opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(0)));
+  prob->addVectorOperator(*opChMold,0, getInvTau());   
+  
+  /// < phi*grad(mu) , grad(psi) >
+  // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
+  Operator *opChLM = new Operator(prob->getFeSpace(0),prob->getFeSpace(1));
+  if (useMobility) {
+    if (doubleWell == 0)
+      opChLM->addSecondOrderTerm(new VecAtQP_SOT(
+	prob->getSolution()->getDOFVector(0),
+	new MobilityCH0(gamma)));
+    else
+      opChLM->addSecondOrderTerm(new VecAtQP_SOT(
+	prob->getSolution()->getDOFVector(0),
+	new MobilityCH1(gamma)));
+  } else
+    opChLM->addSecondOrderTerm(new Simple_SOT(gamma));
+  prob->addMatrixOperator(*opChLM, 0, 1);
+  
+  /// < u_old * grad(c) - u_old * grad(c_old), psi >
+  Operator *uGradC = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+  if (useConservationForm)
+    uGradC->addTerm(new WorldVec_FOT(vel, -1.0), GRD_PSI);
+  else
+    uGradC->addTerm(new WorldVec_FOT(vel), GRD_PHI);
+  uGradC->setUhOld(prob->getSolution()->getDOFVector(0));
+  prob->addVectorOperator(*uGradC, 0);
+  prob->addMatrixOperator(*uGradC, 0, 0);
+  
+  /// < u * grad(c_old) , psi >
+  for (size_t i = 0; i < dow; i++) {
+    Operator *uGradC2 = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
+    if (useConservationForm)
+      uGradC2->addTerm(new VecAndPartialDerivative_FOT(prob->getSolution()->getDOFVector(0), i, -1.0), GRD_PSI);
+    else {
+      uGradC2->addTerm(new PartialDerivative_ZOT(prob->getSolution()->getDOFVector(0), i));
+    }
+    prob->addMatrixOperator(*uGradC2, 0, 2+i);
+  }
+  
+
+  // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
+  // ----------------------------------------------------------------------
+  /// < -3*phi*c*c_old^2 , psi >
+  // -3*c_old^2 * c
+  Operator *opChMPowImpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+    prob->getSolution()->getDOFVector(0),
+    new Pow2Functor(-3.0)));
+  if (doubleWell == 0) {
+    opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
+      prob->getSolution()->getDOFVector(0),
+      new AMDiS::Factor<double>(3.0)));
+    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(-0.5));
+  } else {
+    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
+  }
+  prob->addMatrixOperator(*opChMPowImpl, 1, 0);  
+  
+  /// < -2*phi*c_old^3 , psi >
+  // -2*c_old^3 + 3/2*c_old^2
+  Operator *opChMPowExpl = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+    prob->getSolution()->getDOFVector(0),
+    new Pow3Functor(-2.0)));
+  if (doubleWell == 0) {
+    opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
+      prob->getSolution()->getDOFVector(0),
+      new Pow2Functor(3.0/2.0)));
+  }
+  prob->addVectorOperator(*opChMPowExpl, 1);            
+  
+  /// < -eps^2*phi*grad(c) , grad(psi) >
+  Operator *opChL = new Operator(prob->getFeSpace(1),prob->getFeSpace(0));
+  opChL->addSecondOrderTerm(new Simple_SOT);
+  prob->addMatrixOperator(*opChL, 1, 0, &minusEpsSqr);  
+  
+  /// < phi*mu , psi >
+  Operator *opChM = new Operator(prob->getFeSpace(0),prob->getFeSpace(0));
+  opChM->addZeroOrderTerm(new Simple_ZOT);
+  prob->addMatrixOperator(*opChM, 1, 1);              
+
+  // -----------------------------------------------------------------------------------
+  for (size_t i = 0; i < dow; ++i) {
+    /// < (1/tau)*u'_i , psi >
+    Operator *opTime = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+    opTime->addTerm(new Simple_ZOT(density));
+    prob->addMatrixOperator(*opTime, 2+i, 2+i, getInvTau(), getInvTau());
+    /// < (1/tau)*u_i^old , psi >
+    Operator *opTimeOld = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+    opTime->addTerm(new Phase_ZOT(prob->getSolution()->getDOFVector(2+i), density));
+    prob->addVectorOperator(*opTimeOld, 2+i, getInvTau(), getInvTau());
+ 
+    /// < u^old*grad(u_i^old) , psi >
+    Operator *opUGradU0 = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+    opUGradU0->addTerm(new WorldVector_FOT(velocity, -density), GRD_PHI);
+    opUGradU0->setUhOld(prob->getSolution()->getDOFVector(2+i));
+    if (nonLinTerm == 0) {
+      prob->addVectorOperator(*opUGradU0, 2+i);
+    }
+
+    if (nonLinTerm == 1) {
+      /// < u'*grad(u_i^old) , psi >
+      for (size_t j = 0; j < dow; ++j) {
+        Operator *opUGradU1 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
+        opUGradU1->addTerm(new PartialDerivative_ZOT(
+          prob->getSolution()->getDOFVector(2+i), j, density));
+        prob->addMatrixOperator(*opUGradU1, 2+i, 2+j);
+      }
+    } else if (nonLinTerm == 2) {
+      /// < u^old*grad(u'_i) , psi >
+      for(size_t j = 0; j < dow; ++j) {
+        Operator *opUGradU2 = new Operator(getFeSpace(2+i),getFeSpace(2+i));
+        opUGradU2->addTerm(new VecAndPartialDerivative_FOT(
+          prob->getSolution()->getDOFVector(2+j), j, density), GRD_PHI);
+        prob->addMatrixOperator(*opUGradU2, 2+i, 2+i);
+      }
+    }
+
+//     for (size_t j = 0; j < dow; ++j) {
+//       Operator *opNull = new Operator(getFeSpace(2+i), getFeSpace(2+j));
+//       opNull->addTerm(new Simple_ZOT(0.0));
+//       prob->addMatrixOperator(*opNull, 2+i, 2+j);
+//     }
+
+    /// Diffusion-Operator (including Stress-Tensor for space-dependent viscosity)
+    addLaplaceTerm(i);
+  
+    /// < p , d_i(psi) >
+    Operator *opGradP = new Operator(getFeSpace(2+i),getFeSpace(2+dow));
+    opGradP->addTerm(new PartialDerivative_FOT(i, -1.0), GRD_PSI);
+    prob->addMatrixOperator(*opGradP, 2+i, 2+dow);
+  
+    /// external force, i.e. gravitational force
+    if (norm(force) > DBL_TOL) {
+      Operator *opForce = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+      opForce->addTerm(new Simple_ZOT(force[i]));
+      prob->addVectorOperator(*opForce, 2+i);
+    }
+    
+    /// forces by Cahn-Hilliard terms
+    Operator *opMuGradC = new Operator(getFeSpace(2+i),getFeSpace(0));
+    opMuGradC->addTerm(new VecAndPartialDerivative_FOT(
+      prob->getSolution()->getDOFVector(1), i), GRD_PHI);
+    prob->addMatrixOperator(*opMuGradC, 2+i, 0, &surfaceTension, &surfaceTension);
+    
+    Operator *opMuGradC2 = new Operator(getFeSpace(2+i),getFeSpace(1));
+    opMuGradC2->addTerm(new PartialDerivative_ZOT(
+      prob->getSolution()->getDOFVector(0), i));
+    prob->addMatrixOperator(*opMuGradC2, 2+i, 1, &surfaceTension, &surfaceTension);
+    opMuGradC2->setUhOld(prob->getSolution()->getDOFVector(1));
+    prob->addVectorOperator(*opMuGradC2, 2+i, &surfaceTension, &surfaceTension);
+    
+  }
+
+  /// div(u) = 0
+  for (size_t i = 0; i < dow; ++i) {
+    /// < d_i(u'_i) , psi >
+    Operator *opDivU = new Operator(getFeSpace(2+dow),getFeSpace(2+i));
+    opDivU->addTerm(new PartialDerivative_FOT(i), GRD_PHI);
+    prob->addMatrixOperator(*opDivU, 2+dow, 2+i);
+  }
+}
+
+
+void CahnHilliardNavierStokes::addLaplaceTerm(int i)
+{ FUNCNAME("CahnHilliardNavierStokes::addLaplaceTerm()");
+
+  /// < alpha*[grad(u)+grad(u)^t] , grad(psi) >
+  if (laplaceType == 1) {
+    for (size_t 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));
+      prob->addMatrixOperator(*opLaplaceUi1, 2+i, 2+j);
+    }
+  }
+  
+  /// < alpha*grad(u'_i) , grad(psi) >
+  Operator *opLaplaceUi = new Operator(getFeSpace(2+i), getFeSpace(2+i));
+    opLaplaceUi->addTerm(new Simple_SOT(viscosity));
+  prob->addMatrixOperator(*opLaplaceUi, 2+i, 2+i);
+}

extensions/base_problems/CahnHilliardNavierStokes.h

+/** \file CahnHilliardNavierStokes.h */
+
+#ifndef CAHN_HILLIARD_H
+#define CAHN_HILLIARD_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "ExtendedProblemStat.h"
+#include "HL_SignedDistTraverse.h"
+
+using namespace AMDiS;
+
+class CahnHilliardNavierStokes : public BaseProblem<ExtendedProblemStat>
+{
+public: // definition of types
+
+  typedef BaseProblem<ExtendedProblemStat> super;
+
+public: // public methods
+
+  CahnHilliardNavierStokes(const std::string &name_);
+  ~CahnHilliardNavierStokes();
+
+  virtual void initData();
+
+  double getEpsilon() { return eps; }
+  int getDoubleWellType() { return doubleWell; }
+
+  DOFVector<WorldVector<double> >* getVelocity()
+  {
+    return velocity;
+  }
+  
+protected: // protected methods
+
+  virtual void fillOperators();
+  virtual void fillBoundaryConditions() {}
+  virtual void addLaplaceTerm(int i);
+
+protected: // protected variables
+
+  HL_SignedDistTraverse *reinit;
+
+  DOFVector<WorldVector<double> >* velocity;
+
+  void calcVelocity()
+  {
+    if (dow == 1)
+      transformDOF(prob->getSolution()->getDOFVector(2), velocity, new AMDiS::Vec1WorldVec<double>);
+    else if (dow == 2)
+      transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), velocity, new AMDiS::Vec2WorldVec<double>);
+    else if (dow == 3)
+      transformDOF(prob->getSolution()->getDOFVector(2), prob->getSolution()->getDOFVector(3), prob->getSolution()->getDOFVector(4), velocity, new AMDiS::Vec3WorldVec<double>);
+  }
+  
+  bool useMobility;
+  bool useConservationForm;
+
+  unsigned dim;
+  
+  int laplaceType;
+  int nonLinTerm;
+  int doubleWell;
+
+  // navierStokes parameters
+  double oldTimestep;
+  double viscosity;
+  double density;
+  double c;
+  
+  // Cahn-Hilliard parameters
+  double gamma;
+  double eps;
+  double minusEps;
+  double epsInv;
+  double minusEpsInv;
+  double epsSqr;
+  double minusEpsSqr;
+  
+  double sigma;		// coupling parameter to calculate the surface tension
+  double surfaceTension;// := sigma/epsilon
+  
+  WorldVector<double> force;
+
+  FileVectorWriter *fileWriter;
+};
+
+
+/** \brief
+ * Abstract function for Cahn-Hilliard mobility
+ */
+class MobilityCH0 : public AbstractFunction<double,double>
+{
+  public:
+    MobilityCH0(double gamma_=1.0) : 
+      AbstractFunction<double,double>(4), 
+      gamma(gamma_),
+      delta(1.e-6) { }
+    double operator()(const double &ch) const 
+    {
+      double phase = std::max(0.0, std::min(1.0, ch));
+      double mobility = 0.25*sqr(phase)*sqr(phase-1.0);
+      return gamma * std::max(mobility, delta);
+    }
+
+  protected:
+    double gamma;
+    double delta;
+};
+
+class MobilityCH1 : public AbstractFunction<double,double>
+{
+  public:
+    MobilityCH1(double gamma_=1.0) :
+      AbstractFunction<double,double>(4),
+      gamma(gamma_),
+      delta(1.e-6) { }
+    double operator()(const double &ch) const
+    {
+      double phase = std::max(-1.0, std::min(1.0, ch));
+      double mobility = 0.25*sqr(sqr(phase)-1.0);
+      return gamma * std::max(mobility, delta);
+    }
+
+  protected:
+    double gamma;
+    double delta;
+};
+
+class ScaleFunctor : public AbstractFunction<double,double>
+{
+  public:
+    ScaleFunctor(double factor_=1.0) :
+      AbstractFunction<double,double>(1),
+      factor(factor_) { }
+    double operator()(const double &ch) const
+    {
+      return factor * ch;
+    }
+
+  protected:
+    double factor;
+};
+
+class Pow2Functor : public AbstractFunction<double,double>
+{
+  public:
+    Pow2Functor(double factor_=1.0) :
+      AbstractFunction<double,double>(2),
+      factor(factor_) { }
+    double operator()(const double &ch) const
+    {
+      return factor * sqr(ch);
+    }
+
+  protected:
+    double factor;
+};
+
+class Pow3Functor : public AbstractFunction<double,double>
+{
+  public:
+    Pow3Functor(double factor_=1.0) :
+      AbstractFunction<double,double>(3),
+      factor(factor_) { }
+    double operator()(const double &ch) const
+    {
+      return factor * sqr(ch) * ch;
+    }
+
+  protected:
+    double factor;
+};
+
+#endif // CAHN_HILLIARD_H

extensions/base_problems/CahnHilliard_RB.cc

+#include "CahnHilliard_RB.h"
+#include "Views.h"
+#include "SignedDistFunctors.h"
+#include "PhaseFieldConvert.h"
+#include "Tools.h"
+
+using namespace AMDiS;
+
+CahnHilliard_RB::CahnHilliard_RB(const std::string &name_) :
+  super(name_),
+  useMobility(false),
+  doubleWell(0),
+  gamma(1.0),
+  eps(0.1),
+  minusEps(-0.1),
+  epsInv(10.0),
+  minusEpsInv(-10.0),
+  epsSqr(0.01),
+  minusEpsSqr(-0.01),
+  reinit(NULL)
+{
+  // parameters for CH
+  Parameters::get(name + "->use mobility", useMobility); // mobility
+  Parameters::get(name + "->gamma", gamma); // mobility
+  Parameters::get(name + "->epsilon", eps); // interface width
+  
+  // type of double well: 0= [0,1], 1= [-1,1]
+  Parameters::get(name + "->double-well type", doubleWell); 
+
+  // transformation of the parameters
+  minusEps = -eps;
+  epsInv = 1.0/eps;
+  minusEpsInv = -epsInv;
+  epsSqr = sqr(eps);
+  minusEpsSqr = -epsSqr;
+}
+
+
+CahnHilliard_RB::~CahnHilliard_RB() 
+{ FUNCNAME("CahnHilliard_RB::~CahnHilliard_RB()");
+
+  if (reinit)
+    delete reinit;
+}
+
+
+void CahnHilliard_RB::initData()
+{ FUNCNAME("CahnHilliard_RB::initData()");
+
+  // create instance redistancing class
+  reinit = new HL_SignedDistTraverse("reinit", getMesh()->getDim());
+  
+  super::initData();
+}
+
+
+void CahnHilliard_RB::solveInitialProblem(AdaptInfo *adaptInfo) 
+{ FUNCNAME("CahnHilliard_RB::solveInitialProblem()");
+
+  Flag initFlag = initDataFromFile(adaptInfo);
+
+  if (!initFlag.isSet(DATA_ADOPTED)) {
+    int initialInterface = 0;
+    Initfile::get(name + "->initial interface", initialInterface);