CahnHilliardNavierStokes.h 5.43 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/** \file CahnHilliardNavierStokes.h */

#ifndef CAHN_HILLIARD_NAVIER_STOKES_H
#define CAHN_HILLIARD_NAVIER_STOKES_H

// coupling structures
#include "CouplingIterationInterface.h"
#include "CouplingTimeInterface.h"
#include "CouplingProblemStat.h"

// structures for local refinement
#include "Refinement.h"
#include "MeshFunction_Level.h"

#include "POperators.h"

using namespace AMDiS;

/**
  * \ingroup Problem 
  *
  * \brief
  */
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
template<typename CH_Type = CahnHilliard,
	 typename NS_Type = NavierStokes_TaylorHood>
Praetorius, Simon's avatar
Praetorius, Simon committed
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
class CahnHilliardNavierStokes : public CouplingIterationInterface,
                                 public CouplingTimeInterface,
                                 public CouplingProblemStat
{
public:

  CahnHilliardNavierStokes(std::string name_, CH_Type *chProb_, NS_Type *nsProb_)
  : CouplingProblemStat(name_),
    chProb(chProb_),
    nsProb(nsProb_),
    refFunction(NULL),
    refinement(NULL),
    sigma(0.0),
    surfaceTension(0.0)
  {
    dow = Global::getGeo(WORLD);

    Parameters::get(name + "->sigma", sigma);
    surfaceTension = sigma/chProb->getEpsilon();
  }


  ~CahnHilliardNavierStokes() {
    if (refFunction != NULL)
      delete refFunction;
    if (refinement != NULL)
      delete refinement;
  }


Praetorius, Simon's avatar
Praetorius, Simon committed
56
57
58
59
60
61
62
63
64
65
66
  /**
   * Add the problems to the iterationInterface, timeInterface and couplingProblemStat.
   * As a consequence all problem can be initialized one after another and in the
   * adaption loop they are solved in rotation.
   * At the end the problems are filled with operators and coupling operators as well as
   * boundary conditions are added.
   *
   * In the adaption loop the problems are solved the same order as they are added to the
   * iterationInterface in this method. This order can be changed manually in the oneIteration
   * method.
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
  void initialize(AdaptInfo *adaptInfo)
  {
    for (size_t i = 0; i < chProb->getNumProblems(); i++)
      addProblem(chProb->getProblem(i));
    for (size_t i = 0; i < nsProb->getNumProblems(); i++)
      addProblem(nsProb->getProblem(i));
    
    addIterationInterface(chProb);
    addIterationInterface(nsProb);

    addTimeInterface(chProb);
    addTimeInterface(nsProb);

    CouplingProblemStat::initialize(INIT_ALL);
    dim = getMesh()->getDim();

    fillCouplingOperators();
    // fillOperators and fillBoundaryConditions for chProb and nsProb
    nsProb->initTimeInterface();
    chProb->initTimeInterface();
Praetorius, Simon's avatar
Praetorius, Simon committed
87
//     fillCouplingBoundaryConditions();
Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
90
  }


Praetorius, Simon's avatar
Praetorius, Simon committed
91
92
93
94
95
  /**
   * In this method the operators responsible for the coupling of all problems are defined
   * and added to the problem. The access to the problems is through the baseBroblem->getProblem()
   * method. All baseProblem provide FiniteElemSpaces and solution vectors.
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
  void fillCouplingOperators()
  { FUNCNAME("CahnHilliardNavierStokes::fillCouplingOperators()");
    MSG("CahnHilliardNavierStokes::fillCouplingOperators()");

    for (size_t i = 0; i < dow; i++) {
      // < nu * d_i(c) , theta >
      Operator *opNuGradC = new Operator(nsProb->getFeSpace(i), chProb->getFeSpace(0));
      opNuGradC->addTerm(new VecAndPartialDerivative_ZOT(chProb->getSolution()->getDOFVector(1), chProb->getSolution()->getDOFVector(0), i));
      nsProb->getProblem(0)->addVectorOperator(opNuGradC, i, &surfaceTension, &surfaceTension);
    }

    // < v * grad(c) , theta >
    Operator *opVGradC = new Operator(chProb->getFeSpace(0), chProb->getFeSpace(0));
    opVGradC->addTerm(new WorldVector_FOT(nsProb->getVelocity(), -1.0), GRD_PSI);
    chProb->getProblem()->addMatrixOperator(opVGradC, 0, 0);
  }
  

Praetorius, Simon's avatar
Praetorius, Simon committed
114
115
116
117
118
119
  /**
   * Solves the initial problems.Before this is done the mesh is refined by the Refinement-class
   * This provides local refinement functions for phaseField-refinement or signed-distance-refinement.
   * See MeshFunction_Level.h for some details.
   **/ 
  void solveInitialProblem(AdaptInfo *adaptInfo)
Praetorius, Simon's avatar
Praetorius, Simon committed
120
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
    // meshFunction for klocal refinement around the interface of the phasefield-DOFVector
    refFunction = new PhaseFieldRefinement(chProb->getMesh());

    if (chProb->getDoubleWellType() == 0) {
      refinement = new RefinementLevelDOF(
	chProb->getFeSpace(),
	refFunction,
	chProb->getSolution()->getDOFVector(0)); // phaseField-DOFVector in [0,1]
    } else {
      refinement = new RefinementLevelDOF(
	chProb->getFeSpace(),
	refFunction,
	new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(0))); // phaseField-DOFVector in [-1,1]
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
135

Praetorius, Simon's avatar
Praetorius, Simon committed
136
    // initial refinement
Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
    refinement->refine(0);

Praetorius, Simon's avatar
Praetorius, Simon committed
139
140
141
    for (int i = 0; i < 3; i++) {
      chProb->solveInitialProblem(adaptInfo); 	// update phaseField-DOFVector
      refinement->refine((i < 4 ? 4 : 10));	// do k steps of local refinement
Praetorius, Simon's avatar
Praetorius, Simon committed
142
143
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
    // solve all initial problems of the problems added to the CouplingTimeInterface
    CouplingTimeInterface::solveInitialProblem(adaptInfo);
Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
148
  }


Praetorius, Simon's avatar
Praetorius, Simon committed
149
150
151
152
153
  /**
   * Called at the end of each timestep.
   * If the phase-field has changed the mesh is updated by the refinement-method
   **/
  void closeTimestep(AdaptInfo *adaptInfo)
Praetorius, Simon's avatar
Praetorius, Simon committed
154
155
156
157
158
159
160
161
  {
    CouplingTimeInterface::closeTimestep(adaptInfo);

    refinement->refine(2);
  }

protected:

Praetorius, Simon's avatar
Praetorius, Simon committed
162
163
  CH_Type *chProb;	// CahnHilliard baseProblem
  NS_Type *nsProb;	// Navier-Stokes baseProblem
Praetorius, Simon's avatar
Praetorius, Simon committed
164
165
166
167

  PhaseFieldRefinement* refFunction;
  RefinementLevelDOF *refinement;

Praetorius, Simon's avatar
Praetorius, Simon committed
168
169
  unsigned dim;		// dimension of the mesh
  unsigned dow;		// dimension of the world
Praetorius, Simon's avatar
Praetorius, Simon committed
170

Praetorius, Simon's avatar
Praetorius, Simon committed
171
172
  double sigma;		// coupling parameter to calculate the surface tension
  double surfaceTension;// := sigma/epsilon
Praetorius, Simon's avatar
Praetorius, Simon committed
173
174
175
};

#endif // CAHN_HILLIARD_NAVIER_STOKES_H