PhaseFieldCrystal.h 3.72 KB
Newer Older
1
2
3
4
5
6
/** \file PhaseFieldCrystal.h */

#ifndef PHASE_FIELD_CRYSTAL_H
#define PHASE_FIELD_CRYSTAL_H

#include "AMDiS.h"
Praetorius, Simon's avatar
Praetorius, Simon committed
7
#include "ExtendedProblemStat.h"
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
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
#include "Helpers.h"
#include "POperators.h"

using namespace AMDiS;

const Flag PFC_MESH_ADOPTED = 1<<2;
const Flag PFC_ADOPTED = 1<<3;

/** Phase-field Crystal problem
 */
class PhaseFieldCrystal : public ProblemIterationInterface,
                          public ProblemInstatBase
{
public:

  PhaseFieldCrystal(const std::string &name_);
  ~PhaseFieldCrystal();

  /// Initialisation of the problem.
  virtual void initialize(Flag initFlag,
                  ProblemStat *adoptProblem = NULL,
                  Flag adoptFlag = INIT_NOTHING);
  virtual void initTimeInterface();

  virtual Flag initDataFromFile(AdaptInfo *adaptInfo);
  virtual void solveInitialProblem(AdaptInfo *adaptInfo);
  virtual void transferInitialSolution(AdaptInfo *adaptInfo);

  /// Implementation of ProblemTimeInterface::initTimestep().
  virtual void initTimestep(AdaptInfo *adaptInfo);

  /// Implementation of ProblemTimeInterface::closeTimestep().
  virtual void closeTimestep(AdaptInfo *adaptInfo);

  virtual void beginIteration(AdaptInfo *adaptInfo);
  virtual Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo= FULL_ITERATION);
  virtual void endIteration(AdaptInfo *adaptInfo);

  /// getting methods
  SystemVector *getSolution()
  {
    return prob->getSolution();
  }

  double *getTempParameter() { return &tempParameter; }

  inline Mesh* getMesh(int comp = 0) 
  {
    return prob->getMesh(comp); 
  }

  inline FiniteElemSpace* getFeSpace(int comp = 0) 
  {
    return prob->getFeSpace(comp); 
  }

  std::string getName()
  {
    return name;
  };

  int getNumProblems() 
  {
    return 1;
  };

  ProblemStat *getProblem(int number= 0)
  {
    if (number<0 || number >= getNumProblems())
      throw(std::runtime_error("problem with given number does not exist"));
    if (number == 0)
      return prob;
  };

  ProblemStat *getProblem(std::string name)
  {
    if (name == "pfc")
      return prob;
    else
      throw(std::runtime_error("problem with given name '" + name + "' does not exist"));
  };

  /// setting methods

  void setAssembleMatrixOnlyOnce_butTimestepChange(int i, int j)
  {
    fixedMatrixTimestep.push_back(std::make_pair(i,j));
  }

  void setNumberOfTimesteps(int nTimesteps_)
  {
    nTimesteps= nTimesteps_;
  }

  void writeFiles(AdaptInfo *adaptInfo, bool force);
  void serialize(std::ostream&) {};
  void deserialize(std::istream&) {};

protected:

  virtual void fillOperators();
  virtual void fillBoundaryConditions();

protected:

  ProblemStat *prob;

  bool useMobility;

  unsigned dim;
  unsigned dow;

  int nTimesteps;
  int oldMeshChangeIdx;

  double oldTimestep;
  double tempParameter;
  double r;
  double rho0;
  double density;
  double two;
  double minus2;

private:

  std::vector<std::pair<int,int> > fixedMatrixTimestep;

};


/** \ingroup MainInstat
 * \brief
 * Abstract function to calculate the pure PFC-Energy
 */
class Energy : public BinaryAbstractFunction<double,double,double>
{
  public:
    Energy() : BinaryAbstractFunction<double,double,double>(4) { }

    double operator()(const double &rho, const double &mu) const {
      return -0.25*sqr(sqr(rho)) + 0.5*rho*mu; }
};


class MobilityPfc : public AbstractFunction<double,double>
{
  public:
    MobilityPfc(double density_ = -0.3, double factor_ = 1.0) : 
      AbstractFunction<double,double>(1),
      density(density_),
      factor(factor_),
      delta(1.e-6) { }

    double operator()(const double &rho) const 
    {
      double mobility= abs(rho - 3.0*density)*factor;
      return std::max(mobility, delta);
    }
  
  protected:
    double density;
    double factor;
    double delta;
};

#endif // PHASE_FIELD_CRYSTAL_H