drivenCavity.cc 4.08 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
#include "AMDiS.h"
#include "CahnHilliardNavierStokes.h"
#include "SignedDistFunctors.h"
#include "PhaseFieldConvert.h"
#include "Refinement.h"
#include "MeshFunction_Level.h"

#include "boost/date_time/posix_time/posix_time.hpp"

using namespace AMDiS;
using namespace boost::posix_time;

struct DrivenCavityBC : AbstractFunction<double, WorldVector<double> >
{
  double operator()(const WorldVector<double> &x) const {
    return std::max(0.0, 1.0 - 4.0 * sqr(x[0] - 0.5));
  }
};

class CHNS_DrivenCavity : public CahnHilliardNavierStokes
{
public:
  CHNS_DrivenCavity(std::string name_) : CahnHilliardNavierStokes(name_) {}

  void solveInitialProblem(AdaptInfo* adaptInfo)
  {
    super::solveInitialProblem(adaptInfo);
    /// horizontale Linie
    double a= 0.0, dir= -1.0;
    double initialEps = eps;
    Initfile::get(name + "->initial epsilon", initialEps);
    Initfile::get(name + "->line->pos", a);
    Initfile::get(name + "->line->direction", dir);
    
    /// create phase-field from signed-dist-function
    if (doubleWell == 0) {
      prob->getSolution()->getDOFVector(0)->interpol(new SignedDistFctToPhaseField(initialEps, new Plane(a, dir)));
    } else {
      prob->getSolution()->getDOFVector(0)->interpol(new Plane(a, dir)); 
      transformDOF(prob->getSolution()->getDOFVector(0),
        new SignedDistToCh(initialEps));
    }
  }
  
  void fillBoundaryConditions()
  { FUNCNAME("NS_DrivenCavity::fillBoundaryConditions()");
    
    DOFVector<double> *zeroDOF = new DOFVector<double>(getFeSpace(0), "zero");
    zeroDOF->set(0.0);
    size_t dow = Global::getGeo(WORLD);

    /// at rigid wall: no-slip boundary condition
    for (size_t i = 0; i < dow; i++)
      getProblem(0)->addDirichletBC(1, 2+i, 2+i, zeroDOF);

    /// at upper wall: prescribed velocity
    getProblem(0)->addDirichletBC(2, 2, 2, new DrivenCavityBC);
    getProblem(0)->addDirichletBC(2, 3, 3, zeroDOF);
  }
};


class RefinementTimeInterface : public CouplingTimeInterface
{
public:

  RefinementTimeInterface(CHNS_DrivenCavity *chnsProb_) :
    chnsProb(chnsProb_),
    refFunction(NULL),
    refinement(NULL)
  {
    addTimeInterface(chnsProb);
  }

  ~RefinementTimeInterface()
  {
    if (refFunction)
      delete refFunction;
    if (refinement)
      delete refinement;
  }

  virtual void initTimeInterface()
  {
    chnsProb->initTimeInterface();
  }

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

    refinement->refine(1);
  }

  /// Set initial condition and perform initial refinement
  virtual void solveInitialProblem(AdaptInfo *adaptInfo)
  {
    refFunction= new PhaseFieldRefinement("mesh->refinement",chnsProb->getMesh());
    refinement= new RefinementLevelDOF(
      chnsProb->getProblem(0)->getFeSpace(0),
      refFunction,
      new PhaseDOFView<double>(chnsProb->getProblem(0)->getSolution()->getDOFVector(0)));

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

    bool initialRefinement = true;
    Parameters::get(chnsProb->getName() + "->initial refinement", initialRefinement);
    
    if (initialRefinement) {
      // refine until interfaces is solved
      for (int i = 0; i < 3; ++i) {
	chnsProb->solveInitialProblem(adaptInfo);
	refinement->refine((i < 4 ? 4 : 10));
      }
    }

    CouplingTimeInterface::solveInitialProblem(adaptInfo);
  }

protected:

  CHNS_DrivenCavity *chnsProb;
  PhaseFieldRefinement* refFunction;
  RefinementLevelDOF *refinement;
};


int main(int argc, char** argv)
{ FUNCNAME("main");

  AMDiS::init(argc, argv);

  CHNS_DrivenCavity chnsProb("chns");
  chnsProb.initialize(INIT_ALL);

  RefinementTimeInterface timeInterface(&chnsProb);

  // Adapt-Infos
  AdaptInfo adaptInfo("adapt", chnsProb.getNumComponents());
  AdaptInstationary adaptInstat("adapt", chnsProb, adaptInfo, timeInterface, adaptInfo);

  ptime start_time = microsec_clock::local_time();
  chnsProb.initTimeInterface();
  int error_code = adaptInstat.adapt(); 
  time_duration td = microsec_clock::local_time()-start_time;

  MSG("elapsed time= %d sec\n", td.total_seconds());

  AMDiS::finalize();

  return error_code;
};