navier_stokes.cc 3.87 KB
 Praetorius, Simon committed Mar 19, 2016 1 2 3 4 5 ``````#include #include #include #include `````` Praetorius, Simon committed May 01, 2016 6 ``````#include `````` Praetorius, Simon committed Mar 19, 2016 7 ``````#include `````` Praetorius, Simon committed May 01, 2016 8 ``````#include `````` Praetorius, Simon committed May 07, 2016 9 ``````#include `````` Praetorius, Simon committed Mar 19, 2016 10 `````` `````` 11 12 ``````#ifdef DOW #undef DOW `````` Praetorius, Simon committed Mar 19, 2016 13 ``````#endif `````` 14 ``````#define DOW AMDIS_DOW `````` Praetorius, Simon committed Mar 19, 2016 15 `````` `````` Praetorius, Simon committed Jun 28, 2016 16 17 18 19 ``````#ifndef STRATEGY #define STRATEGY 1 #endif `````` Praetorius, Simon committed Mar 19, 2016 20 21 ``````using namespace AMDiS; `````` Praetorius, Simon committed Jun 28, 2016 22 23 24 ``````using Mesh = Dune::YaspGrid; // using Mesh = Dune::AlbertaGrid; `````` Praetorius, Simon committed Mar 19, 2016 25 ``````// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1 `````` Praetorius, Simon committed Nov 19, 2017 26 ``````using StokesParam = LagrangeTraits; `````` Praetorius, Simon committed Mar 19, 2016 27 ``````using StokesProblem = ProblemStat; `````` Praetorius, Simon committed May 01, 2016 28 ``````using StokesProblemInstat = ProblemInstat; `````` Praetorius, Simon committed Mar 19, 2016 29 30 31 32 `````` int main(int argc, char** argv) { AMDiS::init(argc, argv); `````` Praetorius, Simon committed Jun 28, 2016 33 `````` `````` Praetorius, Simon committed Mar 19, 2016 34 35 `````` StokesProblem prob("stokes"); prob.initialize(INIT_ALL); `````` Praetorius, Simon committed Jun 28, 2016 36 `````` `````` Praetorius, Simon committed May 01, 2016 37 38 `````` StokesProblemInstat probInstat("stokes", prob); probInstat.initialize(INIT_UH_OLD); `````` Praetorius, Simon committed Jun 28, 2016 39 `````` `````` Praetorius, Simon committed Mar 19, 2016 40 41 42 43 44 45 `````` double viscosity = 1.0; double density = 1.0; double vel = 1.0; Parameters::get("stokes->viscosity", viscosity); Parameters::get("stokes->density", density); Parameters::get("stokes->boundary velocity", vel); `````` Praetorius, Simon committed Jun 28, 2016 46 `````` `````` 47 `````` AdaptInfo adaptInfo("adapt", prob.getNumComponents()); `````` Praetorius, Simon committed Jun 28, 2016 48 `````` `````` Praetorius, Simon committed Mar 19, 2016 49 `````` // define the differential operators `````` Praetorius, Simon committed Nov 11, 2017 50 `````` using Op = StokesProblem::ElementOperator; `````` Praetorius, Simon committed Nov 09, 2017 51 `````` forEach(range_<0,DOW>, [&](auto const _i) `````` Praetorius, Simon committed Jun 28, 2016 52 `````` { `````` Praetorius, Simon committed May 01, 2016 53 `````` // <1/tau * u_i, v_i> `````` Praetorius, Simon committed Jun 28, 2016 54 55 `````` auto opTime = Op::zot( density ); auto opTimeOld = Op::zot( valueOf(prob.getSolution(_i), density) ); `````` Praetorius, Simon committed Nov 10, 2017 56 57 `````` prob.addMatrixOperator(*opTime, _i, _i, probInstat.getInvTau()); prob.addVectorOperator(*opTimeOld, _i, probInstat.getInvTau()); `````` Praetorius, Simon committed Jun 28, 2016 58 59 `````` #if STRATEGY == 1 `````` Praetorius, Simon committed May 01, 2016 60 `````` // <(u * nabla)u_i^old, v_i> `````` Praetorius, Simon committed Nov 09, 2017 61 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed May 01, 2016 62 `````` { `````` Praetorius, Simon committed Jun 28, 2016 63 `````` auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) ); `````` Praetorius, Simon committed Nov 10, 2017 64 `````` prob.addMatrixOperator(*opNonlin, _i, _j); `````` Praetorius, Simon committed May 01, 2016 65 `````` }); `````` Praetorius, Simon committed Jun 28, 2016 66 ``````#elif STRATEGY == 2 `````` Praetorius, Simon committed May 01, 2016 67 `````` // <(u^old * nabla)u_i, v_i> `````` Praetorius, Simon committed Nov 09, 2017 68 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed May 01, 2016 69 `````` { `````` Praetorius, Simon committed Jun 28, 2016 70 `````` auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI ); `````` Praetorius, Simon committed Nov 10, 2017 71 `````` prob.addMatrixOperator(*opNonlin, _i, _i); `````` Praetorius, Simon committed Jun 28, 2016 72 73 74 `````` }); #else // <(u^old * nabla)u_i^old, v_i> `````` Praetorius, Simon committed Nov 09, 2017 75 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed Jun 28, 2016 76 77 `````` { auto opNonlin = Op::zot( density * valueOf(prob.getSolution(_j)) * derivativeOf(prob.getSolution(_i), _j) ); `````` Praetorius, Simon committed Nov 10, 2017 78 `````` prob.addVectorOperator(*opNonlin, _i); `````` Praetorius, Simon committed May 01, 2016 79 `````` }); `````` Praetorius, Simon committed Jun 28, 2016 80 81 82 83 84 85 86 87 88 89 90 ``````#endif auto opL = std::make_pair( Op::sot( viscosity ), false ); // auto opP = std::make_pair( Op::fot( 1.0, _i, GRD_PSI ), false); // auto opDiv = std::make_pair( Op::fot( 1.0, _i, GRD_PHI ), false); // const int i = _i; prob.addMatrixOperator({{ {{i, i}, opL}, {{i,DOW}, opP}, {{DOW,i}, opDiv} }}); `````` Praetorius, Simon committed Mar 19, 2016 91 `````` }); `````` Praetorius, Simon committed Jun 28, 2016 92 93 94 95 96 `````` auto opZero = std::make_pair( Op::zot( 0.0 ), false ); prob.addMatrixOperator({{ {{DOW, DOW}, opZero} }}); `````` Praetorius, Simon committed Mar 19, 2016 97 98 99 `````` // define boundary regions auto left = [](auto const& x) { return x[0] < 1.e-8; }; auto box = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; }; `````` Praetorius, Simon committed Jun 28, 2016 100 101 `````` auto corner = [](auto const& x) { return x[0] == 0.0 && x[1] == 0.0; }; `````` Praetorius, Simon committed Mar 19, 2016 102 103 104 `````` // define boundary values auto zero = [](auto const& x) { return 0.0; }; auto parabolic = [vel](auto const& x) { return vel*4.0*x[1]*(1.0 - x[1]); }; `````` Praetorius, Simon committed Jun 28, 2016 105 `````` `````` Praetorius, Simon committed Mar 19, 2016 106 107 108 `````` // set boundary conditions for (size_t i = 0; i < DOW; ++i) prob.addDirichletBC(box, i, i, zero); `````` Praetorius, Simon committed Jun 28, 2016 109 `````` `````` Praetorius, Simon committed Mar 19, 2016 110 111 `````` prob.addDirichletBC(left, 0, 0, zero); prob.addDirichletBC(left, 1, 1, parabolic); `````` Praetorius, Simon committed Jun 28, 2016 112 113 114 115 `````` // pressure fixation prob.addDirichletBC(corner, DOW, DOW, zero); `````` Praetorius, Simon committed Mar 20, 2016 116 `````` *prob.getSolution() = 0.0; `````` Praetorius, Simon committed Jun 28, 2016 117 `````` `````` Praetorius, Simon committed May 01, 2016 118 119 `````` AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo); adapt.adapt(); `````` Praetorius, Simon committed Jun 28, 2016 120 `````` `````` Praetorius, Simon committed Mar 19, 2016 121 `````` // output solution `````` Praetorius, Simon committed Apr 27, 2016 122 `````` prob.writeFiles(adaptInfo); `````` Praetorius, Simon committed Jun 28, 2016 123 `````` `````` Praetorius, Simon committed Mar 19, 2016 124 125 126 `````` AMDiS::finalize(); return 0; }``````