CahnHilliard_.cc 6.93 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
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
#include "CahnHilliard_.h"
// #include "Views.h"
#include "SignedDistFunctors.h"
#include "PhaseFieldConvert.h"

using namespace AMDiS;

CahnHilliard_::CahnHilliard_(const std::string &name_, bool createProblem) :
  super(name_, createProblem),
  useMobility(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); 

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


CahnHilliard_::~CahnHilliard_() 
{ FUNCNAME("CahnHilliard_::~CahnHilliard_()"); }


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(1)->interpol(new Plane(a, dir));
    }
    else if (initialInterface == 1) {
      /// schraege Linie
      double theta = m_pi/4.0;
      prob->getSolution()->getDOFVector(1)->interpol(new PlaneRotation(0.0, theta, 1.0));
      transformDOFInterpolation(prob->getSolution()->getDOFVector(1),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(1)->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(1)->interpol(new Plane(a, -1.0));
      transformDOFInterpolation(prob->getSolution()->getDOFVector(1),new Plane(b, 1.0), new AMDiS::Max<double>);
    }
    else if (initialInterface == 4) {
      /// Kreis
      double radius= 1.0, x=0, y=0;
      Initfile::get(name + "->circle->radius", radius);
      Initfile::get(name + "->circle->center_x", x);
      Initfile::get(name + "->circle->center_y", y);
      WorldVector<double> w;
      w[0]=x; w[1]=y+adaptInfo->getTime();
      prob->getSolution()->getDOFVector(1)->interpol(new Circle(radius, w));
    } 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(1)->interpol(new Rectangle(width, height, center));
    }

    /// create phase-field from signed-dist-function
    if (doubleWell == 0) {
118
119
      forEachDOF(prob->getSolution()->getDOFVector(1),
        new SignedDistToPhaseField(initialEps, 1.0/(2.0*sqrt(2.0))));
120
    } else {
121
122
      forEachDOF(prob->getSolution()->getDOFVector(1),
        new SignedDistToCh(initialEps, 1.0/sqrt(2.0)));
123
124
125
126
127
128
129
130
131
132
    }
  }
}


void CahnHilliard_::fillOperators()
{ FUNCNAME("CahnHilliard_::fillOperators()");
  MSG("CahnHilliard_::fillOperators()\n");

  const FiniteElemSpace* feSpace = prob->getFeSpace(0);
133
  int degree = feSpace->getBasisFcts()->getDegree();
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
  
  // c
  Operator *opChMnew = new Operator(feSpace,feSpace);
  opChMnew->addZeroOrderTerm(new Simple_ZOT);
  Operator *opChMold = new Operator(feSpace,feSpace);
  opChMold->addZeroOrderTerm(new VecAtQP_ZOT(prob->getSolution()->getDOFVector(1)));
  // -nabla*(grad(c))
  Operator *opChL = new Operator(feSpace,feSpace);
  opChL->addSecondOrderTerm(new Simple_SOT);
  
  // div(M(c)grad(mu)), with M(c)=gamma/4*(c^2-1)^2
  Operator *opChLM = new Operator(feSpace,feSpace);
  opChLM->addSecondOrderTerm(new Simple_SOT(gamma));
  
  // -2*c_old^3 + 3/2*c_old^2
  Operator *opChMPowExpl = new Operator(feSpace,feSpace);
  opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
    prob->getSolution()->getDOFVector(1),
152
    new AMDiS::Pow<3>(-2.0,3*degree)));
153
154
155
  if (doubleWell == 0) {
    opChMPowExpl->addZeroOrderTerm(new VecAtQP_ZOT(
      prob->getSolution()->getDOFVector(1),
156
      new AMDiS::Pow<2>(3.0/2.0,2*degree)));
157
158
159
160
161
162
  }

  // -3*c_old^2 * c
  Operator *opChMPowImpl = new Operator(feSpace,feSpace);
  opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
    prob->getSolution()->getDOFVector(1),
163
    new AMDiS::Pow<2>(-3.0,2*degree)));
164
165
166
  if (doubleWell == 0) {
    opChMPowImpl->addZeroOrderTerm(new VecAtQP_ZOT(
      prob->getSolution()->getDOFVector(1),
167
      new AMDiS::Factor<double>(3.0,degree)));
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
    opChMPowImpl->addZeroOrderTerm(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]
  // ----------------------------------------------------------------------
  prob->addMatrixOperator(*opChMnew,0,0);              /// < phi*mu , psi >
  
  prob->addMatrixOperator(*opChMPowImpl,0,1);          /// < -3*phi*c*c_old^2 , psi >
  prob->addMatrixOperator(*opChL,0,1, &minusEpsSqr);   /// < -eps^2*phi*grad(c) , grad(psi) >
  // . . . vectorOperators . . . . . . . . . . . . . . .
  prob->addVectorOperator(*opChMPowExpl,0);            /// < -2*phi*c_old^3 , psi >

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