CahnHilliard.cc 8.02 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
/******************************************************************************
 *
 * 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.
 * 
 ******************************************************************************/
17
#include "CahnHilliard.h"
18
// #include "Views.h"
19
20
21
#include "SignedDistFunctors.h"
#include "PhaseFieldConvert.h"

Praetorius, Simon's avatar
Praetorius, Simon committed
22
23
24
#include "HL_SignedDistTraverse.h"
#include "Recovery.h"

25
26
using namespace AMDiS;

Praetorius, Simon's avatar
Praetorius, Simon committed
27
28
29
30
namespace details {

template<typename P>
CahnHilliard<P>::CahnHilliard(const std::string &name_) :
31
32
  super(name_),
  useMobility(false),
Praetorius, Simon's avatar
Praetorius, Simon committed
33
  useReinit(false),
34
35
36
37
38
39
40
  doubleWell(0),
  gamma(1.0),
  eps(0.1),
  minusEps(-0.1),
  epsInv(10.0),
  minusEpsInv(-10.0),
  epsSqr(0.01),
41
  minusEpsSqr(-0.01)
42
43
{
  // parameters for CH
44
45
46
  Parameters::get(name_ + "->use mobility", useMobility); // mobility
  Parameters::get(name_ + "->gamma", gamma); // mobility
  Parameters::get(name_ + "->epsilon", eps); // interface width
47
48
  
  // type of double well: 0= [0,1], 1= [-1,1]
49
  Parameters::get(name_ + "->double-well type", doubleWell); 
Praetorius, Simon's avatar
Praetorius, Simon committed
50
  
Praetorius, Simon's avatar
Praetorius, Simon committed
51
  Parameters::get(name_ + "->use reinit", useReinit);
52
53
54
55
56
57
58
59
60
61

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


Praetorius, Simon's avatar
Praetorius, Simon committed
62
63
64
65
template<typename P>
void CahnHilliard<P>::solveInitialProblem(AdaptInfo *adaptInfo) 
{
  Flag initFlag = self::initDataFromFile(adaptInfo);
66
67
68

  if (!initFlag.isSet(DATA_ADOPTED)) {
    int initialInterface = 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
69
    Initfile::get(self::name + "->initial interface", initialInterface);
70
    double initialEps = eps;
Praetorius, Simon's avatar
Praetorius, Simon committed
71
    Initfile::get(self::name + "->initial epsilon", initialEps);
72
73
74
75

    if (initialInterface == 0) {
      /// horizontale Linie
      double a= 0.0, dir= -1.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
78
      Initfile::get(self::name + "->line->pos", a);
      Initfile::get(self::name + "->line->direction", dir);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
79
80
81
82
    }
    else if (initialInterface == 1) {
      /// schraege Linie
      double theta = m_pi/4.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
83
84
      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>);
85
86
87
88
    }
    else if (initialInterface == 2) {
      /// Ellipse
      double a= 1.0, b= 1.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
89
90
91
      Initfile::get(self::name + "->ellipse->a", a);
      Initfile::get(self::name + "->ellipse->b", b);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Ellipse(a,b));
92
93
94
95
    }
    else if (initialInterface == 3) {
      /// zwei horizontale Linien
      double a= 0.75, b= 0.375;
Praetorius, Simon's avatar
Praetorius, Simon committed
96
97
98
99
      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>);
100
101
102
103
    }
    else if (initialInterface == 4) {
      /// Kreis
      double radius= 1.0;
Praetorius, Simon's avatar
Praetorius, Simon committed
104
105
      Initfile::get(self::name + "->kreis->radius", radius);
      self::prob->getSolution()->getDOFVector(0)->interpol(new Circle(radius));
106
107
108
109
110
    } else if (initialInterface == 5) {
      /// Rechteck
      double width = 0.5;
      double height = 0.3;
      WorldVector<double> center; center.set(0.5);
Praetorius, Simon's avatar
Praetorius, Simon committed
111
112
113
114
      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));
115
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
116
117
118
119
    
    // TODO: Redistancing einfügen!
    if (useReinit) {
      FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace(
Praetorius, Simon's avatar
Praetorius, Simon committed
120
121
122
	const_cast<DOFAdmin*>(self::getMesh()->getVertexAdmin()),
	Lagrange::getLagrange(self::getMesh()->getDim(), 1),
	self::getMesh(),
Praetorius, Simon's avatar
Praetorius, Simon committed
123
124
	"P1");
      DOFVector<double> tmp(feSpace, "tmp");
Praetorius, Simon's avatar
Praetorius, Simon committed
125
      tmp.interpol(self::prob->getSolution()->getDOFVector(0));
Praetorius, Simon's avatar
Praetorius, Simon committed
126
    
Praetorius, Simon's avatar
Praetorius, Simon committed
127
      HL_SignedDistTraverse reinit("reinit", self::getMesh()->getDim());
Praetorius, Simon's avatar
Praetorius, Simon committed
128
129
      reinit.calcSignedDistFct(adaptInfo, &tmp);
    
130
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
Praetorius, Simon's avatar
Praetorius, Simon committed
131
      Recovery recovery(L2_NORM, 1);
Praetorius, Simon's avatar
Praetorius, Simon committed
132
      recovery.recoveryUh(&tmp, *self::prob->getSolution()->getDOFVector(0));
133
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
134
      self::prob->getSolution()->getDOFVector(0)->interpol(&tmp);
135
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
136
137
    }
    
138
139
140

    /// create phase-field from signed-dist-function
    if (doubleWell == 0) {
Praetorius, Simon's avatar
Praetorius, Simon committed
141
      forEachDOF(self::prob->getSolution()->getDOFVector(0),
142
143
        new SignedDistToPhaseField(initialEps));
    } else {
Praetorius, Simon's avatar
Praetorius, Simon committed
144
      forEachDOF(self::prob->getSolution()->getDOFVector(0),
145
146
147
148
149
150
        new SignedDistToCh(initialEps));
    }
  }
}


Praetorius, Simon's avatar
Praetorius, Simon committed
151
152
153
154
template<typename P>
void CahnHilliard<P>::fillOperators()
{  
  const FiniteElemSpace* feSpace = self::prob->self::getFeSpace();
155
  int degree = feSpace->getBasisFcts()->getDegree();
156
157

  // c
158
  Operator *opChMnew = new Operator(feSpace, feSpace);
159
  opChMnew->addTerm(new Simple_ZOT);
160
  Operator *opChMold = new Operator(feSpace, feSpace);
Praetorius, Simon's avatar
Praetorius, Simon committed
161
  opChMold->addTerm(new VecAtQP_ZOT(self::prob->getSolution()->getDOFVector(0)));
162
  // -nabla*(grad(c))
163
  Operator *opChL = new Operator(feSpace, feSpace);
164
  opChL->addTerm(new Simple_SOT);
165
166
  
  // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
167
  Operator *opChLM = new Operator(feSpace, feSpace);
168
169
  if (useMobility) {
    if (doubleWell == 0)
170
      opChLM->addTerm(new VecAtQP_SOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
171
	self::prob->getSolution()->getDOFVector(0),
172
	new MobilityCH0(gamma, degree)));
173
    else
174
      opChLM->addTerm(new VecAtQP_SOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
175
	self::prob->getSolution()->getDOFVector(0),
176
	new MobilityCH1(gamma, degree)));
177
  } else
178
    opChLM->addTerm(new Simple_SOT(gamma));
179
180
  
  // -2*c_old^3 + 3/2*c_old^2
181
  Operator *opChMPowExpl = new Operator(feSpace, feSpace);
182
  opChMPowExpl->addTerm(new VecAtQP_ZOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
183
    self::prob->getSolution()->getDOFVector(0),
184
    new AMDiS::Pow<3>(-2.0, 3*degree)));
185
  if (doubleWell == 0) {
186
    opChMPowExpl->addTerm(new VecAtQP_ZOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
187
      self::prob->getSolution()->getDOFVector(0),
188
      new AMDiS::Pow<2>(3.0/2.0, 2*degree)));
189
190
191
  }

  // -3*c_old^2 * c
192
  Operator *opChMPowImpl = new Operator(feSpace, feSpace);
193
  opChMPowImpl->addTerm(new VecAtQP_ZOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
194
    self::prob->getSolution()->getDOFVector(0),
195
    new AMDiS::Pow<2>(-3.0, 2*degree)));
196
  if (doubleWell == 0) {
197
    opChMPowImpl->addTerm(new VecAtQP_ZOT(
Praetorius, Simon's avatar
Praetorius, Simon committed
198
      self::prob->getSolution()->getDOFVector(0),
199
      NULL, 3.0));
200
    opChMPowImpl->addTerm(new Simple_ZOT(-0.5));
201
202
203
204
205
206
  } else {
    opChMPowImpl->addZeroOrderTerm(new Simple_ZOT(1.0));
  }

  // mu + eps^2*laplace(c) + c - 3*(c_old^2)*c = -2*c_old^3 [+ BC]
  // ----------------------------------------------------------------------
Praetorius, Simon's avatar
Praetorius, Simon committed
207
208
209
  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 >
210
  // . . . vectorOperators . . . . . . . . . . . . . . .
Praetorius, Simon's avatar
Praetorius, Simon committed
211
  self::prob->addVectorOperator(*opChMPowExpl,0);            /// < -2*phi*c_old^3 , psi >
212

213
214
215
  
  // dt(c) = laplace(mu) - u*grad(c)
  // -----------------------------------
Praetorius, Simon's avatar
Praetorius, Simon committed
216
217
  self::prob->addMatrixOperator(*opChMnew,1,0, self::getInvTau()); /// < phi*c/tau , psi >
  self::prob->addMatrixOperator(*opChLM,1,1);                /// < phi*grad(mu) , grad(psi) >
218
  // . . . vectorOperators . . . . . . . . . . . . . . .
Praetorius, Simon's avatar
Praetorius, Simon committed
219
  self::prob->addVectorOperator(*opChMold,1, self::getInvTau());   /// < phi*c^old/tau , psi >
220
  
221
}
222
223


Praetorius, Simon's avatar
Praetorius, Simon committed
224
225
226
227
228
template<typename P>
void CahnHilliard<P>::finalizeData()
{  
  self::setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
  self::setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
229
  if (!useMobility)
Praetorius, Simon's avatar
Praetorius, Simon committed
230
    self::setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
231
232
}

Praetorius, Simon's avatar
Praetorius, Simon committed
233
} // end namespace details