CahnHilliardNavierStokes.h 4.56 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
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
/** \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"
#include "CouplingOperators.h"

using namespace AMDiS;

/**
  * \ingroup Problem 
  *
  * \brief
  */
template<typename CH_Type=CahnHilliard, typename NS_Type=NavierStokes_TaylorHood>
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;
  }


  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();
    fillCouplingBoundaryConditions();
  }


  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);
// 
//       // stabilizing term
//       Operator *stabilBeltrami = new Operator(nsProb->getFeSpace(i),  nsProb->getFeSpace(i));
//       stabilBeltrami->addTerm(new MatrixGradient_SOT(
//         chProb->getSolution()->getDOFVector(0), new ProjectionMatrix(surfaceTension), new DivFct())); // factor = sigma ?
//       nsProb->getProblem(0)->addMatrixOperator(stabilBeltrami, i, i, nsProb->getTau(), nsProb->getTau());
    }

    // < 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);
  }


  void fillCouplingBoundaryConditions()
  { FUNCNAME("CahnHilliardNavierStokes::fillCouplingBoundaryConditions()");
  
  }
  

  /// Solves the initial problem.
  virtual void solveInitialProblem(AdaptInfo *adaptInfo) 
  {
    refFunction= new PhaseFieldRefinement(chProb->getMesh());
    refinement= new RefinementLevelDOF(
      chProb->getFeSpace(),
      refFunction,
      new PhaseDOFView<double>(chProb->getSolution()->getDOFVector(0))); // phaseField-DOFVector

    // set initial values
    refinement->refine(0);

    for (int i= 0; i < 3; ++i) {
      chProb->solveInitialProblem(adaptInfo); // initial phaseField
      refinement->refine((i < 4 ? 4 : 10));
    }
    
    CouplingTimeInterface::solveInitialProblem(adaptInfo);
  }


  /// Called at the beginning of each timestep
  virtual void initTimestep(AdaptInfo *adaptInfo) 
  {
    CouplingTimeInterface::initTimestep(adaptInfo);
  }


  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
  {
    CouplingIterationInterface::oneIteration(adaptInfo, toDo);
  }


  /// Called at the end of each timestep.
  virtual void closeTimestep(AdaptInfo *adaptInfo) 
  {
    CouplingTimeInterface::closeTimestep(adaptInfo);

    refinement->refine(2);
  }

protected:

  CH_Type *chProb;
  NS_Type *nsProb;

  PhaseFieldRefinement* refFunction;
  RefinementLevelDOF *refinement;

  unsigned dim;
  unsigned dow;

  double sigma;
  double surfaceTension;
};

#endif // CAHN_HILLIARD_NAVIER_STOKES_H