demo for polarization ordering

extensions/base_problems/PolarizationField.h

+/******************************************************************************
+ *
+ * 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.
+ * 
+ ******************************************************************************/
+
+#ifndef POLARIZATION_FIELD_H
+#define POLARIZATION_FIELD_H
+
+#include "AMDiS.h"
+#include "BaseProblem.h"
+#include "ExtendedProblemStat.h"
+#include "GenericOperatorTerm.h"
+
+
+namespace detail 
+{
+  using namespace AMDiS;
+  
+  /** \ingroup PolarizationField
+  * \brief
+  * Simulation of the relaxation of an orientation field (polarization field)
+  */
+  
+  template<typename ProblemStatType>
+  class PolarizationField : public BaseProblem<ProblemStatType>
+  {
+  public: // typedefs
+    
+    typedef PolarizationField<ProblemStatType> self;
+    typedef BaseProblem<ProblemStatType> super;
+    
+  public: // methods
+
+    PolarizationField(const std::string &name_);
+
+    ~PolarizationField();
+
+    void initData() override;
+
+    void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+
+    void closeTimestep(AdaptInfo *adaptInfo) override;
+    
+    void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
+
+    // === getting/setting methods ===
+
+    DOFVector<WorldVector<double> >* getVectorField()
+    { FUNCNAME_DBG("getVectorField()");
+      TEST_EXIT_DBG(vectorField != NULL)
+	("vectorField is NULL!");
+      return vectorField;
+    }
+    
+    DOFVector<double> *getOldSolution(int i)
+    { FUNCNAME_DBG("getOldSolution()");
+      TEST_EXIT_DBG(oldSolution[i] != NULL)
+	("Index out of range, or oldSolution not provided for this component!\n");
+      return oldSolution[i];
+    }
+
+    void fillOperators() override;
+    virtual void fillLaplacian();
+
+  protected: // variables
+
+    DOFVector<WorldVector<double> >* vectorField;
+
+    void calcVectorField()
+    {
+      if (self::dow == 1)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    vectorField, new AMDiS::Vec1WorldVec<double>);
+      else if (self::dow == 2)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    self::prob->getSolution()->getDOFVector(1), 
+		    vectorField, new AMDiS::Vec2WorldVec<double>);
+      else if (self::dow == 3)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    self::prob->getSolution()->getDOFVector(1), 
+		    self::prob->getSolution()->getDOFVector(2), 
+		    vectorField, new AMDiS::Vec3WorldVec<double>);
+    }
+
+    double oldTimestep;
+    double minus1;
+    
+    double alpha2;
+    double alpha4;
+    double epsilon;
+    double epsInv;
+    double K;
+
+    FileVectorWriter *fileWriter;
+    std::vector<DOFVector<double>*> oldSolution;
+  };
+
+} // end namespace detail
+
+#include "PolarizationField.hh"
+
+typedef ::detail::PolarizationField<AMDiS::ProblemStat> PolarizationField;
+
+
+
+#endif // POLARIZATION_FIELD_H
+

extensions/base_problems/PolarizationField.hh

+/******************************************************************************
+ *
+ * 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 "Helpers.h"
+#include "POperators.h"
+
+namespace detail {
+
+using namespace AMDiS;
+
+template<typename P> 
+PolarizationField<P>::PolarizationField(const std::string &name_) :
+  super(name_, true),
+  vectorField(NULL),
+  oldTimestep(0.0),
+  minus1(-1.0),
+  alpha2(1.0),
+  alpha4(1.0),
+  epsilon(0.1),
+  K(1.0),
+  fileWriter(NULL)
+{
+  oldSolution.resize(self::dow);
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i] = NULL;
+  
+  Parameters::get(self::name + "->alpha2", alpha2);
+  Parameters::get(self::name + "->alpha4", alpha4);
+  
+  Parameters::get(self::name + "->epsilon", epsilon);
+  epsInv = 1.0/epsilon;
+  
+  Parameters::get(self::name + "->K", K);
+}
+
+template<typename P> 
+PolarizationField<P>::~PolarizationField() 
+{   
+  if (vectorField != NULL) {
+    delete vectorField;
+    vectorField = NULL;
+  }
+  
+  for (size_t i = 0; i < self::dow; i++) {
+    if (oldSolution[i] != NULL)
+      delete oldSolution[i];
+    oldSolution[i] = NULL;
+  }
+ 
+  if (fileWriter) { 
+    delete fileWriter;
+    fileWriter = NULL;
+  }
+}
+
+
+template<typename P> 
+void PolarizationField<P>::initData()
+{ 
+  if (vectorField == NULL)
+    vectorField = new DOFVector<WorldVector<double> >(self::getFeSpace(0), "vectorField");
+  
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i] = new DOFVector<double>(self::getFeSpace(i), "old(v_"+Helpers::toString(i)+")");
+  
+  fileWriter = new FileVectorWriter(self::name + "->vectorField->output", self::getFeSpace()->getMesh(), vectorField);
+
+  super::initData();
+}
+
+
+template<typename P> 
+void PolarizationField<P>::transferInitialSolution(AdaptInfo *adaptInfo)
+{ 
+  calcVectorField();
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i]->copy(*self::prob->getSolution()->getDOFVector(i));
+  
+  super::transferInitialSolution(adaptInfo);
+}
+
+
+template<typename P> 
+void PolarizationField<P>::fillOperators()
+{ 
+  WorldVector<DOFVector<double>* > vec;
+  for (size_t k = 0; k < self::dow; k++)
+    vec[k] = self::prob->getSolution()->getDOFVector(k);
+  
+  const FiniteElemSpace* feSpace = self::getFeSpace(0);
+  
+  // fill operators for component P
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < (1/tau)*P_i , psi >
+    Operator *opTime = new Operator(feSpace, feSpace);
+    addZOT(opTime, constant(1.0));
+    self::prob->addMatrixOperator(*opTime, i, i, self::getInvTau(), self::getInvTau());
+    
+    /// < (1/tau)*P_i^old , psi >
+    Operator *opTimeOld = new Operator(feSpace, feSpace);
+    addZOT(opTimeOld, valueOf(getOldSolution(i)));
+    self::prob->addVectorOperator(*opTimeOld, i, self::getInvTau(), self::getInvTau());
+
+    /// Diffusion-Operator
+    Operator *opLaplace = new Operator(feSpace, feSpace);
+    addSOT(opLaplace, constant(alpha2));
+    addZOT(opLaplace, constant(alpha4));
+    self::prob->addMatrixOperator(*opLaplace, i, i+self::dow);
+  }
+
+  // fill operators for component P#
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < P# , psi >
+    Operator *opM = new Operator(feSpace, feSpace);
+    addZOT(opM, constant(1.0));
+    self::prob->addMatrixOperator(*opM, i+self::dow, i+self::dow);
+    
+    /// < (-C1 - C4*P^2)*P , psi >
+    Operator *opNonlin = new Operator(feSpace, feSpace);
+    addZOT(opNonlin, epsInv * (1.0 - pow<2>(valueOf(vectorField))) );
+    self::prob->addMatrixOperator(*opNonlin, i+self::dow, i);
+    
+    for (size_t j = 0; j < self::dow; ++j) {
+      Operator *opNonlin2 = new Operator(feSpace, feSpace);
+      addZOT(opNonlin2, (-2.0*epsInv) * valueOf(vec[i]) * valueOf(vec[j]));
+      self::prob->addMatrixOperator(*opNonlin, i+self::dow, j);
+    }
+    
+    Operator *opNonlin3 = new Operator(feSpace, feSpace);
+    addZOT(opNonlin3, (-2.0*epsInv) * valueOf(vec[i]) * pow<2>(valueOf(vectorField)) );
+    self::prob->addVectorOperator(*opNonlin3, i+self::dow);
+  }
+  
+  fillLaplacian();
+
+  Operator *opNull = new Operator(feSpace, feSpace);
+  addZOT(opNull, constant(0.0));
+  
+  for (size_t i = 0; i < self::dow; ++i) {
+    for (size_t j = i+1; j < self::dow; ++j) {
+      self::prob->addMatrixOperator(*opNull, i, j);
+      self::prob->addMatrixOperator(*opNull, j, i);
+      self::prob->addMatrixOperator(*opNull, i, j+self::dow);
+      self::prob->addMatrixOperator(*opNull, j, i+self::dow);
+      self::prob->addMatrixOperator(*opNull, i+self::dow, j+self::dow);
+      self::prob->addMatrixOperator(*opNull, j+self::dow, i+self::dow);
+    }
+  }
+}
+
+
+template<typename P> 
+void PolarizationField<P>::fillLaplacian()
+{
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < -K*grad(P) , grad(psi) >
+    Operator *opL = new Operator(feSpace, feSpace);
+    addSOT(opL, constant(-K));
+    self::prob->addMatrixOperator(*opL, i+self::dow, i);
+  }
+}
+
+
+
+template<typename P> 
+void PolarizationField<P>::closeTimestep(AdaptInfo *adaptInfo)
+{ FUNCNAME("PolarizationField::closeTimestep()");
+
+  calcVectorField();
+  for (size_t i = 0; i < self::dow; i++)
+    oldSolution[i]->copy(*self::prob->getSolution()->getDOFVector(i));
+  
+  super::closeTimestep(adaptInfo);
+}
+
+
+template<typename P> 
+void PolarizationField<P>::writeFiles(AdaptInfo *adaptInfo, bool force)
+{ FUNCNAME("PolarizationField::closeTimestep()");
+
+  super::writeFiles(adaptInfo, force);
+  self::fileWriter->writeFiles(adaptInfo, false);
+}
+
+} // end namespace detail
\ No newline at end of file

extensions/base_problems/PolarizationField_RB.h

+/******************************************************************************
+ *
+ * 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.
+ * 
+ ******************************************************************************/
+
+#ifndef POLARIZATION_FIELD_RB_H
+#define POLARIZATION_FIELD_RB_H
+
+#include "AMDiS.h"
+#include "BaseProblem_RB.h"
+#include "ExtendedProblemStat.h"
+#include "GenericOperatorTerm.h"
+
+
+namespace detail 
+{
+  using namespace AMDiS;
+  
+  /** \ingroup PolarizationField_RB
+  * \brief
+  * Simulation of the relaxation of an orientation field (polarization field)
+  */
+  
+  template<typename ProblemStatType>
+  class PolarizationField_RB : public BaseProblem_RB //<ProblemStatType>
+  {
+  public: // typedefs
+    
+    typedef PolarizationField_RB<ProblemStatType> self;
+    typedef BaseProblem_RB/*<ProblemStatType>*/ super;
+    
+  public: // methods
+
+    PolarizationField_RB(const std::string &name_);
+
+    ~PolarizationField_RB();
+
+    void initData() override;
+
+    void transferInitialSolution(AdaptInfo *adaptInfo) override;  
+
+    void closeTimestep(AdaptInfo *adaptInfo) override;
+    
+    void writeFiles(AdaptInfo *adaptInfo, bool force = false) override;
+
+    // === getting/setting methods ===
+
+    DOFVector<WorldVector<double> >* getVectorField()
+    { FUNCNAME_DBG("getVectorField()");
+      TEST_EXIT_DBG(vectorField != NULL)
+	("vectorField is NULL!");
+      return vectorField;
+    }
+
+    void fillOperators() override;
+    virtual void fillLaplacian();
+
+  protected: // variables
+
+    DOFVector<WorldVector<double> >* vectorField;
+
+    void calcVectorField()
+    {
+      if (self::dow == 1)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    vectorField, new AMDiS::Vec1WorldVec<double>);
+      else if (self::dow == 2)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    self::prob->getSolution()->getDOFVector(1), 
+		    vectorField, new AMDiS::Vec2WorldVec<double>);
+      else if (self::dow == 3)
+	transformDOF(self::prob->getSolution()->getDOFVector(0), 
+		    self::prob->getSolution()->getDOFVector(1), 
+		    self::prob->getSolution()->getDOFVector(2), 
+		    vectorField, new AMDiS::Vec3WorldVec<double>);
+    }
+
+    double minus1;
+    
+    double alpha2;
+    double alpha4;
+    double epsilon;
+    double epsInv;
+    double K;
+
+    FileVectorWriter *fileWriter;
+  };
+
+} // end namespace detail
+
+#include "PolarizationField_RB.hh"
+
+typedef ::detail::PolarizationField_RB<AMDiS::ProblemStat> PolarizationField_RB;
+
+
+
+#endif // POLARIZATION_FIELD_RB_H
+

extensions/base_problems/PolarizationField_RB.hh

+/******************************************************************************
+ *
+ * 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 "Helpers.h"
+#include "POperators.h"
+
+namespace detail {
+
+using namespace AMDiS;
+
+template<typename P> 
+PolarizationField_RB<P>::PolarizationField_RB(const std::string &name_) :
+  super(name_, true),
+  vectorField(NULL),
+  minus1(-1.0),
+  alpha2(1.0),
+  alpha4(1.0),
+  epsilon(0.1),
+  K(1.0),
+  fileWriter(NULL)
+{  
+  Parameters::get(self::name + "->alpha2", alpha2);
+  Parameters::get(self::name + "->alpha4", alpha4);
+  
+  Parameters::get(self::name + "->epsilon", epsilon);
+  epsInv = 1.0/epsilon;
+  
+  Parameters::get(self::name + "->K", K);
+}
+
+template<typename P> 
+PolarizationField_RB<P>::~PolarizationField_RB() 
+{   
+  if (vectorField != NULL) {
+    delete vectorField;
+    vectorField = NULL;
+  }
+ 
+  if (fileWriter) { 
+    delete fileWriter;
+    fileWriter = NULL;
+  }
+}
+
+
+template<typename P> 
+void PolarizationField_RB<P>::initData()
+{ 
+  if (vectorField == NULL)
+    vectorField = new DOFVector<WorldVector<double> >(self::getFeSpace(0), "vectorField");
+  
+  fileWriter = new FileVectorWriter(self::name + "->vectorField->output", self::getFeSpace()->getMesh(), vectorField);
+
+  super::initData();
+}
+
+
+template<typename P> 
+void PolarizationField_RB<P>::transferInitialSolution(AdaptInfo *adaptInfo)
+{ 
+  calcVectorField();  
+  super::transferInitialSolution(adaptInfo);
+}
+
+
+template<typename P> 
+void PolarizationField_RB<P>::fillOperators()
+{ 
+  WorldVector<DOFVector<double>* > vec;
+  WorldVector<DOFVector<double>* > vec0;
+  for (size_t k = 0; k < self::dow; k++) {
+    vec[k] = self::prob->getStageSolution(k);
+    vec0[k] = self::prob->getUnVec(k);
+  }
+  
+  const FiniteElemSpace* feSpace = self::getFeSpace(0);
+  
+  // fill operators for component P
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < (1/tau)*P_i , psi >
+    self::prob->addTimeOperator(i, i);	// dt(P_i)
+
+    /// Diffusion-Operator
+    Operator *opLaplace = new Operator(feSpace, feSpace);
+    opLaplace->setUhOld(self::prob->getStageSolution(i+self::dow)); // P#_i
+    addSOT(opLaplace, constant(-alpha2));
+    addZOT(opLaplace, constant(-alpha4));
+    
+    self::prob->addMatrixOperator(*opLaplace, i, i+self::dow, &minus1, &minus1);
+    self::prob->addVectorOperator(*opLaplace, i);
+  }
+
+  // fill operators for component P#
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < P# , psi >
+    Operator *opM = new Operator(feSpace, feSpace);
+    opM->setUhOld(self::prob->getStageSolution(i+self::dow)); // P#_i
+    addZOT(opM, constant(1.0));
+    
+    self::prob->addMatrixOperator(*opM, i+self::dow, i+self::dow, &minus1, &minus1);
+    self::prob->addVectorOperator(*opM, i+self::dow);
+    
+    /// < (-C1 - C4*P^2)*P , psi >
+    Operator *opNonlinExpl = new Operator(feSpace, feSpace);
+    Operator *opNonlinImpl = new Operator(feSpace, feSpace);
+    if (self::dow == 2) {
+      addZOT(opNonlinExpl, epsInv * (1.0 - pow<2>(valueOf(vec[0])) - pow<2>(valueOf(vec[1]))) );
+      addZOT(opNonlinImpl, epsInv * (1.0 - pow<2>(valueOf(vec0[0])) - pow<2>(valueOf(vec0[1]))) ); // (1 - |P|^2)*dP
+    } else if (self::dow == 3) {
+      addZOT(opNonlinExpl, epsInv * (1.0 - pow<2>(valueOf(vec[0])) - pow<2>(valueOf(vec[1])) - pow<2>(valueOf(vec[2]))) );
+      addZOT(opNonlinImpl, epsInv * (1.0 - pow<2>(valueOf(vec0[0])) - pow<2>(valueOf(vec0[1])) - pow<2>(valueOf(vec0[2]))) );
+    }
+    self::prob->addMatrixOperator(*opNonlinImpl, i+self::dow, i, &minus1, &minus1);
+    self::prob->addVectorOperator(*opNonlinExpl, i+self::dow);
+        
+    for (size_t j = 0; j < self::dow; ++j) {
+      Operator *opNonlinImpl2 = new Operator(feSpace, feSpace);
+      addZOT(opNonlinImpl2, (-2.0*epsInv) * valueOf(vec0[i]) * valueOf(vec0[j]));
+      self::prob->addMatrixOperator(*opNonlinImpl2, i+self::dow, j, &minus1, &minus1);
+    }
+  }
+
+  Operator *opNull = new Operator(feSpace, feSpace);
+  addZOT(opNull, constant(0.0));
+  
+  for (size_t i = 0; i < self::dow; ++i) {
+    for (size_t j = i+1; j < self::dow; ++j) {
+      self::prob->addMatrixOperator(*opNull, i, j);
+      self::prob->addMatrixOperator(*opNull, j, i);
+      self::prob->addMatrixOperator(*opNull, i, j+self::dow);
+      self::prob->addMatrixOperator(*opNull, j, i+self::dow);
+      self::prob->addMatrixOperator(*opNull, i+self::dow, j+self::dow);
+      self::prob->addMatrixOperator(*opNull, j+self::dow, i+self::dow);
+    }
+  }
+}
+
+
+template<typename P> 
+void PolarizationField_RB<P>::fillLaplacian()
+{
+  for (size_t i = 0; i < self::dow; ++i) {
+    /// < -K*grad(P) , grad(psi) >
+    Operator *opLaplace2 = new Operator(feSpace, feSpace);
+    opLaplace2->setUhOld(self::prob->getStageSolution(i)); // P#_i
+    addSOT(opLaplace2, constant(-K));
+    
+    self::prob->addMatrixOperator(*opLaplace2, i+self::dow, i, &minus1, &minus1);