stokes1.cc 2.69 KB
Newer Older
1 2 3 4 5 6
#include <iostream>
#include <ctime>
#include <cmath>

#include <dune/amdis/AMDiS.hpp>
#include <dune/amdis/ProblemStat.hpp>
7
#include <dune/amdis/Operators.hpp>
8 9 10 11 12 13 14 15 16 17

#ifdef DOW
#undef DOW
#endif
#define DOW AMDIS_DOW

using namespace AMDiS;

// 3 components: velocity with polynomial degree 2 and pressure with polynomial degree 1

18
using Grid = Dune::YaspGrid<AMDIS_DIM>;
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
using StokesParam   = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;

int main(int argc, char** argv)
{
    AMDiS::init(argc, argv);

    StokesProblem prob("stokes");
    prob.initialize(INIT_ALL);

    double viscosity = 1.0;
    Parameters::get("stokes->viscosity", viscosity);

    // tree-paths for components
    auto _v = 0_c;
    auto _p = 1_c;

36
    // <viscosity*grad(u_i), grad(v_i)>
37 38
    for (std::size_t i = 0; i < DOW; ++i) {
      auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
39 40 41
      prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));
    }

42 43 44
    // <d_i(v_i), p>
    auto opP = makeOperator(tag::divtestvec_trial{}, 1.0);
    prob.addMatrixOperator(opP, _v, _p);
45

46 47
    // <q, d_i(u_i)>
    auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0);
48 49
    prob.addMatrixOperator(opDiv, _p, _v);

50
    auto opZero = makeOperator(tag::test_trial{}, 0.0);
51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
    prob.addMatrixOperator(opZero, _p, _p);


    // define boundary regions
    auto left = [](auto const& x) { return x[0] < 1.e-8; };
    auto not_left = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };

    // define boundary values
    auto parabolic_y = [](auto const& x) -> Dune::FieldVector<double,DOW>
    {
      return {0.0, x[1]*(1.0 - x[1])};
    };

    auto zero = [](auto const& x) -> Dune::FieldVector<double,DOW>
    {
      return {0.0, 0.0};
    };

    // set boundary conditions for velocity
    prob.addDirichletBC(left, _v, _v, parabolic_y);
    prob.addDirichletBC(not_left, _v, _v, zero);

73
    prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, _p, _p, 0.0);
74 75 76 77 78

    AdaptInfo adaptInfo("adapt", prob.getNumComponents());

    // assemble and solve system
    prob.buildAfterCoarsen(adaptInfo, Flag(0));
79 80 81 82 83 84 85 86 87 88 89

    // write matrix to file
    if (Parameters::get<int>("stokesMesh->global refinements").value_or(0) < 2) {
      mtl::io::matrix_market_ostream out("stokes_matrix.mtx");
      out << prob.getSystemMatrix()->getMatrix();

      std::cout << "A = \n" << prob.getSystemMatrix()->getMatrix() << '\n';
      std::cout << "x = \n" << prob.getSolution()->getVector() << '\n';
      std::cout << "b = \n" << prob.getRhs()->getVector() << '\n';
    }

90 91 92 93 94 95 96 97
    prob.solve(adaptInfo);

    // output solution
    prob.writeFiles(adaptInfo);

    AMDiS::finalize();
    return 0;
}