navierStokes.cc 3.05 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
#include "AMDiS.h"
#include "NavierStokes_TaylorHood_RB.h"
#include "navierStokes.h"
#include "time/ExtendedRosenbrockAdaptInstationary.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;



class NS_Channel : public NavierStokes_TaylorHood_RB
{
public:
  NS_Channel(std::string name_) : NavierStokes_TaylorHood_RB(name_) {}

  void solveInitialProblem(AdaptInfo *adaptInfo)
  {
    super::solveInitialProblem(adaptInfo);

    size_t nVertices = 0;
    WorldVector<double> x;
    Parameters::get("obstacle->num vertices",nVertices);
    std::vector<WorldVector<double> > v(nVertices,x);
    for (size_t i = 0; i < nVertices; i++)
      Parameters::get("obstacle->vertex["+boost::lexical_cast<std::string>(i)+"]",v[i]);
    v.push_back(v[0]);

    SignedDistRefinement refFunction(getMesh());
    RefinementLevelCoords2 refinement(
      getFeSpace(),
      &refFunction,
      new Polygon(v));

    // initial refinement
    refinement.refine(10);
  }
protected:
  void fillBoundaryConditions()
  { FUNCNAME("NS_DrivenCavity::fillBoundaryConditions()");
    
    AbstractFunction<double, WorldVector<double> > *zero = new AMDiS::Const<double, WorldVector<double> >(0.0);
    size_t dow = Global::getGeo(WORLD);

    // +------ 5 ------+
    // |               |
    // 2   # <--1      3
    // |               |
    // +------ 4 ------+


//     x[0] = 2.0;
//     x[1] = 2.0;
#if DEBUG != 0
  DOFVector<double>* polygon = new DOFVector<double>(getFeSpace(0), "polygon");
  polygon->interpol(new Polygon(v));
  VtkWriter::writeFile(polygon, "polygon.vtu");
#endif
    
    /// at rigid wall: no-slip boundary condition
    for (size_t i = 0; i < dow; i++) {
//       getProblem(0)->addImplicitDirichletBC(*(new Polygon(v)), i, i, *zero);
//       getProblem(0)->addSingularDirichletBC(x, i, i, *zero);
      getProblem(0)->addDirichletBC(1, i, i, zero);
      getProblem(0)->addDirichletBC(4, i, i, zero);
      getProblem(0)->addDirichletBC(5, i, i, zero);
    }

    double H = 4.1;
    double Um = 1.5;
    Parameters::get("mesh->H",H);
    Parameters::get("ns->Um",Um);

    /// at left wall: prescribed velocity
    getProblem(0)->addDirichletBC(2, 0, 0, new InflowBC(H,Um));
    getProblem(0)->addDirichletBC(2, 1, 1, zero);
  }
};


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

  AMDiS::init(argc, argv);

  NS_Channel nsProb("ns");
  nsProb.initialize(INIT_ALL | INIT_EXACT_SOLUTION);

  // Adapt-Infos
  AdaptInfo adaptInfo("adapt", nsProb.getNumComponents());

  // adaption loop - solve ch-prob and ns-prob
//   AdaptInstationary adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);
  ExtendedRosenbrockAdaptInstationary<NS_Channel> adaptInstat("adapt", nsProb, adaptInfo, nsProb, adaptInfo);

  ptime start_time = microsec_clock::local_time();
  nsProb.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;
};