CahnHilliard.cc 7.72 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
27
28
29
using namespace AMDiS;

CahnHilliard::CahnHilliard(const std::string &name_) :
  super(name_),
  useMobility(false),
Praetorius, Simon's avatar
Praetorius, Simon committed
30
  useReinit(false),
31
32
33
34
35
36
37
  doubleWell(0),
  gamma(1.0),
  eps(0.1),
  minusEps(-0.1),
  epsInv(10.0),
  minusEpsInv(-10.0),
  epsSqr(0.01),
38
  minusEpsSqr(-0.01)
39
40
{
  // parameters for CH
41
42
43
  Parameters::get(name_ + "->use mobility", useMobility); // mobility
  Parameters::get(name_ + "->gamma", gamma); // mobility
  Parameters::get(name_ + "->epsilon", eps); // interface width
44
45
  
  // type of double well: 0= [0,1], 1= [-1,1]
46
  Parameters::get(name_ + "->double-well type", doubleWell); 
Praetorius, Simon's avatar
Praetorius, Simon committed
47
48
  
  Parameters::get(name + "->use reinit", useReinit);
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

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


void CahnHilliard::solveInitialProblem(AdaptInfo *adaptInfo) 
{ FUNCNAME("CahnHilliard::solveInitialProblem()");

  Flag initFlag = initDataFromFile(adaptInfo);

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

    if (initialInterface == 0) {
      /// horizontale Linie
      double a= 0.0, dir= -1.0;
      Initfile::get(name + "->line->pos", a);
      Initfile::get(name + "->line->direction", dir);
      prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir));
    }
    else if (initialInterface == 1) {
      /// schraege Linie
      double theta = m_pi/4.0;
      prob->getSolution()->getDOFVector(0)->interpol(new PlaneRotation(0.0, theta, 1.0));
      transformDOFInterpolation(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(name + "->ellipse->a", a);
      Initfile::get(name + "->ellipse->b", b);
      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(name + "->lines->pos1", a);
      Initfile::get(name + "->lines->pos2", b);
      prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, -1.0));
      transformDOFInterpolation(prob->getSolution()->getDOFVector(0),new Plane(b, 1.0), new AMDiS::Max<double>);
    }
    else if (initialInterface == 4) {
      /// Kreis
      double radius= 1.0;
      Initfile::get(name + "->kreis->radius", radius);
      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(name + "->rectangle->width", width);
      Initfile::get(name + "->rectangle->height", height);
      Initfile::get(name + "->rectangle->center", center);
      prob->getSolution()->getDOFVector(0)->interpol(new Rectangle(width, height, center));
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    
    // TODO: Redistancing einfügen!
    if (useReinit) {
      FiniteElemSpace* feSpace = FiniteElemSpace::provideFeSpace(
	const_cast<DOFAdmin*>(getMesh()->getVertexAdmin()),
	Lagrange::getLagrange(getMesh()->getDim(), 1),
	getMesh(),
	"P1");
      DOFVector<double> tmp(feSpace, "tmp");
      tmp.interpol(prob->getSolution()->getDOFVector(0));
    
      HL_SignedDistTraverse reinit("reinit", getMesh()->getDim());
      reinit.calcSignedDistFct(adaptInfo, &tmp);
    
127
#ifndef HAVE_PARALLEL_DOMAIN_AMDIS
Praetorius, Simon's avatar
Praetorius, Simon committed
128
129
      Recovery recovery(L2_NORM, 1);
      recovery.recoveryUh(&tmp, *prob->getSolution()->getDOFVector(0));
130
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
131
      prob->getSolution()->getDOFVector(0)->interpol(&tmp);
132
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
133
134
    }
    
135
136
137

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


void CahnHilliard::fillOperators()
{ FUNCNAME("CahnHilliard::fillOperators()");
150
  
151
152
  const FiniteElemSpace* feSpace = prob->getFeSpace();
  int degree = feSpace->getBasisFcts()->getDegree();
153
154

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

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

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

210
211
212
213
214
215
216
217
  
  // dt(c) = laplace(mu) - u*grad(c)
  // -----------------------------------
  prob->addMatrixOperator(*opChMnew,1,0, getInvTau()); /// < phi*c/tau , psi >
  prob->addMatrixOperator(*opChLM,1,1);                /// < phi*grad(mu) , grad(psi) >
  // . . . vectorOperators . . . . . . . . . . . . . . .
  prob->addVectorOperator(*opChMold,1, getInvTau());   /// < phi*c^old/tau , psi >
  
218
}
219
220
221
222
223
224
225
226
227
228
229


void CahnHilliard::finalizeData()
{ FUNCNAME("CahnHilliard::finalizeData()");
  
  setAssembleMatrixOnlyOnce_butTimestepChange(0,1);
  setAssembleMatrixOnlyOnce_butTimestepChange(1,0);
  if (!useMobility)
    setAssembleMatrixOnlyOnce_butTimestepChange(1,1);
}