navier_stokes.cc 3.94 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 Nov 15, 2017 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 15, 2017 56 57 `````` prob.addMatrixOperator(opTime, _i, _i, probInstat.getInvTau()); prob.addVectorOperator(opTimeOld, _i, probInstat.getInvTau()); `````` Praetorius, Simon committed Jun 28, 2016 58 `````` `````` Praetorius, Simon committed May 01, 2016 59 `````` // <(u * nabla)u_i^old, v_i> `````` Praetorius, Simon committed Nov 09, 2017 60 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed May 01, 2016 61 `````` { `````` Praetorius, Simon committed Jun 28, 2016 62 `````` auto opNonlin = Op::zot( derivativeOf(prob.getSolution(_i), _j, density) ); `````` Praetorius, Simon committed Nov 15, 2017 63 `````` prob.addMatrixOperator(opNonlin, _i, _j); `````` Praetorius, Simon committed May 01, 2016 64 `````` }); `````` Praetorius, Simon committed Nov 15, 2017 65 `````` `````` Praetorius, Simon committed May 01, 2016 66 `````` // <(u^old * nabla)u_i, v_i> `````` Praetorius, Simon committed Nov 09, 2017 67 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed May 01, 2016 68 `````` { `````` Praetorius, Simon committed Jun 28, 2016 69 `````` auto opNonlin = Op::fot( valueOf(prob.getSolution(_j), density), _j, GRD_PHI ); `````` Praetorius, Simon committed Nov 15, 2017 70 `````` prob.addMatrixOperator(opNonlin, _i, _i); `````` Praetorius, Simon committed Jun 28, 2016 71 `````` }); `````` Praetorius, Simon committed Nov 15, 2017 72 `````` `````` Praetorius, Simon committed Jun 28, 2016 73 `````` // <(u^old * nabla)u_i^old, v_i> `````` Praetorius, Simon committed Nov 09, 2017 74 `````` forEach(range_<0, DOW>, [&](auto const _j) `````` Praetorius, Simon committed Jun 28, 2016 75 `````` { `````` Praetorius, Simon committed Nov 15, 2017 76 77 78 79 80 `````` auto opNonlin = Op::zot( func([density](double uj, double dj_ui) { return density * uj * dj_ui; }, valueOf(prob.getSolution(_j)), derivativeOf(prob.getSolution(_i), _j)) ); prob.addVectorOperator(opNonlin, _i); `````` Praetorius, Simon committed May 01, 2016 81 `````` }); `````` Praetorius, Simon committed Jun 28, 2016 82 `````` `````` Praetorius, Simon committed Nov 15, 2017 83 84 85 `````` 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); // `````` Praetorius, Simon committed Jun 28, 2016 86 87 88 89 90 91 `````` const int i = _i; prob.addMatrixOperator({{ {{i, i}, opL}, {{i,DOW}, opP}, {{DOW,i}, opDiv} }}); `````` Praetorius, Simon committed Mar 19, 2016 92 `````` }); `````` Praetorius, Simon committed Jun 28, 2016 93 94 95 96 97 `````` auto opZero = std::make_pair( Op::zot( 0.0 ), false ); prob.addMatrixOperator({{ {{DOW, DOW}, opZero} }}); `````` Praetorius, Simon committed Mar 19, 2016 98 99 100 `````` // 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 101 102 `````` auto corner = [](auto const& x) { return x[0] == 0.0 && x[1] == 0.0; }; `````` Praetorius, Simon committed Mar 19, 2016 103 104 105 `````` // 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 106 `````` `````` Praetorius, Simon committed Mar 19, 2016 107 108 `````` // set boundary conditions for (size_t i = 0; i < DOW; ++i) `````` Praetorius, Simon committed Nov 15, 2017 109 `````` prob.addDirichletBC(box, i, i, zero); `````` Praetorius, Simon committed Jun 28, 2016 110 `````` `````` Praetorius, Simon committed Mar 19, 2016 111 112 `````` prob.addDirichletBC(left, 0, 0, zero); prob.addDirichletBC(left, 1, 1, parabolic); `````` Praetorius, Simon committed Jun 28, 2016 113 114 115 116 `````` // pressure fixation prob.addDirichletBC(corner, DOW, DOW, zero); `````` Praetorius, Simon committed Mar 20, 2016 117 `````` *prob.getSolution() = 0.0; `````` Praetorius, Simon committed Jun 28, 2016 118 `````` `````` Praetorius, Simon committed May 01, 2016 119 120 `````` AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo); adapt.adapt(); `````` Praetorius, Simon committed Jun 28, 2016 121 `````` `````` Praetorius, Simon committed Mar 19, 2016 122 `````` // output solution `````` Praetorius, Simon committed Apr 27, 2016 123 `````` prob.writeFiles(adaptInfo); `````` Praetorius, Simon committed Jun 28, 2016 124 `````` `````` Praetorius, Simon committed Mar 19, 2016 125 126 127 `````` AMDiS::finalize(); return 0; }``````