CahnHilliard.cc

/******************************************************************************
 *
 * Extension of AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: Simon Praetorius et al.
 *
 * 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.
 * 
 ******************************************************************************/
#include "CahnHilliard.h"
// #include "Views.h"
#include "SignedDistFunctors.h"
#include "PhaseFieldConvert.h"

#include "HL_SignedDistTraverse.h"
#include "Recovery.h"

using namespace AMDiS;

namespace details {

template<typename P>
CahnHilliard<P>::CahnHilliard(const std::string &name_) :
  super(name_),
  useMobility(false),
  useReinit(false),
  doubleWell(0),
  gamma(1.0),
  eps(0.1),
  minusEps(-0.1),
  epsInv(10.0),
  minusEpsInv(-10.0),
  epsSqr(0.01),
  minusEpsSqr(-0.01)
{
  // 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::get(name_ + "->use reinit", useReinit);

  // transformation of the parameters
  minusEps = -eps;
  epsInv = 1.0/eps;
  minusEpsInv = -epsInv;
  epsSqr = sqr(eps);
  minusEpsSqr = -epsSqr;
}


template<typename P>
void CahnHilliard<P>::solveInitialProblem(AdaptInfo *adaptInfo) 
{
  Flag initFlag = self::initDataFromFile(adaptInfo);

  if (!initFlag.isSet(DATA_ADOPTED)) {
    int initialInterface = 0;
    Initfile::get(self::name + "->initial interface", initialInterface);
    double initialEps = eps;
    Initfile::get(self::name + "->initial epsilon", initialEps);

    if (initialInterface == 0) {
      /// horizontale Linie
      double a= 0.0, dir= -1.0;
      Initfile::get(self::name + "->line->pos", a);
      Initfile::get(self::name + "->line->direction", dir);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
    }
    else if (initialInterface == 1) {
      /// schraege Linie
      double theta = m_pi/4.0;
      self::prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new PlaneRotation(0.0, -theta, -1.0), new AMDiS::Min<double>);
    }
    else if (initialInterface == 2) {
      /// Ellipse
      double a= 1.0, b= 1.0;
      Initfile::get(self::name + "->ellipse->a", a);
      Initfile::get(self::name + "->ellipse->b", b);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b));
    }
    else if (initialInterface == 3) {
      /// zwei horizontale Linien
      double a= 0.75, b= 0.375;
      Initfile::get(self::name + "->lines->pos1", a);
      Initfile::get(self::name + "->lines->pos2", b);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
      transformDOFInterpolation(self::prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
    }
    else if (initialInterface == 4) {
      /// Kreis
      double radius= 1.0;
      Initfile::get(self::name + "->kreis->radius", radius);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius));
    } else if (initialInterface == 5) {
      /// Rechteck
      double width = 0.5;
      double height = 0.3;
      WorldVector<double> center; center.set(0.5);
      Initfile::get(self::name + "->rectangle->width", width);
      Initfile::get(self::name + "->rectangle->height", height);
      Initfile::get(self::name + "->rectangle->center", center);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center));
    }
    
    // TODO: Redistancing einfügen!
    if (useReinit) {
      FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace(
	const_cast<DOFAdmin*>(self::getMesh()->getVertexAdmin()),
	Lagrange::getLagrange(self::getMesh()->getDim(), 1),
	self::getMesh(),
	"P1");
      DOFVector<double> tmp(feSpace, "tmp");
      tmp.interpol(self::prob->getSolution()->getDOFVector(0));
    
      HL_SignedDistTraverse reinit("reinit", self::getMesh()->getDim());
      reinit.calcSignedDistFct(adaptInfo, &tmp);
    
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
      Recovery recovery(L2_NORM, 1);
      recovery.recoveryUh(&tmp, *self::prob->getSolution()->getDOFVector(0));
#else
      self::prob->getSolution()->getDOFVector(0)->interpol(&tmp);
#endif
    }
    

    /// create phase-field from signed-dist-function
    if (doubleWell == 0) {
      forEachDOF(self::prob->getSolution()->getDOFVector(0),
        new SignedDistToPhaseField(initialEps));
    } else {
      forEachDOF(self::prob->getSolution()->getDOFVector(0),
        new SignedDistToCh(initialEps));
    }
  }
}


template<typename P>
void CahnHilliard<P>::fillOperators()
{  
  const FiniteElemSpace* feSpace = self::prob->self::getFeSpace();
  int degree = feSpace->getBasisFcts()->getDegree();

  // c
  Operator *opChMnew = new Operator(feSpace, feSpace);
  opChMnew->addTerm(new Simple_ZOT);
  Operator *opChMold = new Operator(feSpace, feSpace);
  opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
  // -nabla*(grad(c))
  Operator *opChL = new Operator(feSpace, feSpace);
  opChL->addTerm(new Simple_SOT);
  
  // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
  Operator *opChLM = new Operator(feSpace, feSpace);
  if (useMobility) {
    if (doubleWell == 0)
      opChLM->addTerm(new VecAtQP_SOT(
	self::prob->getSolution()->getDOFVector(0),
	new MobilityCH0(gamma, degree)));
    else
      opChLM->addTerm(new VecAtQP_SOT(
	self::prob->getSolution()->getDOFVector(0),
	new MobilityCH1(gamma, degree)));
  } else
    opChLM->addTerm(new Simple_SOT(gamma));
  
  // -2*c_old^3 + 3/2*c_old^2
  Operator *opChMPowExpl = new Operator(feSpace, feSpace);
  opChMPowExpl->addTerm(new VecAtQP_ZOT(
    self::prob->getSolution()->getDOFVector(0),
    new AMDiS::Pow<3>(-2.0, 3*degree)));
  if (doubleWell == 0) {
    opChMPowExpl->addTerm(new VecAtQP_ZOT(
      self::prob->getSolution()->getDOFVector(0),
      new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
  }

  // -3*c_old^2 * c
  Operator *opChMPowImpl = new Operator(feSpace, feSpace);
  opChMPowImpl->addTerm(new VecAtQP_ZOT(
    self::prob->getSolution()->getDOFVector(0),
    new AMDiS::Pow<2>(-3.0, 2*degree)));
  if (doubleWell == 0) {
    opChMPowImpl->addTerm(new VecAtQP_ZOT(
      self::prob->getSolution()->getDOFVector(0),
      NULL, 3.0));
    opChMPowImpl->addTerm(new Simple_ZOT(-0.5));
  } else {
    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
  }

  // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
  // ----------------------------------------------------------------------
  self::prob->addMatrixOperator(*opChMPowImpl,0,0);          /// < -3*phi*c*c_old^2 , psi >
  self::prob->addMatrixOperator(*opChL,0,0, &minusEpsSqr);   /// < -eps^2*phi*grad(c) , grad(psi) >
  self::prob->addMatrixOperator(*opChMnew,0,1);              /// < phi*mu , psi >
  // . . . vectorOperators . . . . . . . . . . . . . . .
  self::prob->addVectorOperator(*opChMPowExpl,0);            /// < -2*phi*c_old^3 , psi >

  
  // dt(c) = laplace(mu) - u*grad(c)
  // -----------------------------------
  self::prob->addMatrixOperator(*opChMnew,1,0, self::getInvTau()); /// < phi*c/tau , psi >
  self::prob->addMatrixOperator(*opChLM,1,1);                /// < phi*grad(mu) , grad(psi) >
  // . . . vectorOperators . . . . . . . . . . . . . . .
  self::prob->addVectorOperator(*opChMold,1, self::getInvTau());   /// < phi*c^old/tau , psi >
  
}


template<typename P>
void CahnHilliard<P>::finalizeData()
{  
  self::setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
  self::setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
  if (!useMobility)
    self::setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
}

} // end namespace details