pfc.cc

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:

#include <iostream>
#include <ctime>
#include <cmath>

#include <amdis/AMDiS.hpp>
#include <amdis/AdaptInstationary.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/ProblemInstat.hpp>
#include <amdis/Terms.hpp>

#include "MTLPfcPrecon.hpp"

using namespace AMDiS;

// using Mesh = Dune::YaspGrid<AMDIS_DIM>;
using Mesh = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;

using PfcParam = LagrangeBasis<Mesh::LeafGridView, 2, 2, 2>;


// 3 component with polynomial degree 1
using PfcProblem = ProblemStat<PfcParam>;
using PfcProblemInstat = ProblemInstat<PfcParam>;

using MatrixType = PfcProblem::SystemMatrixType::MultiMatrix;
using VectorType = PfcProblem::SystemVectorType::BaseVector;

using Precon = MTLPfcPrecon<MatrixType,VectorType>;

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

    auto creator = new Precon::Creator;
    CreatorMap<typename Precon::PreconBase>::addCreator("pfc", creator);

    PfcProblem prob("pfc");
    prob.initialize(INIT_ALL);

    PfcProblemInstat probInstat("pfc", prob);
    probInstat.initialize(INIT_UH_OLD);

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

    double r = -0.4;
    double psi_mean = -0.3;
    double M0 = 1.0;

    auto& Mu  = prob.getSolution<0>();
    auto& Psi = prob.getSolution<1>();
    auto& Nu  = prob.getSolution<2>();

    using Op = PfcProblem::ElementOperator;

    // < [-(1+r) - 3*psi^2]*u, v > + < 2*grad(u), grad(v) >
    auto opLhs01_ = Op::zot( -(1.0 + r) );
    opLhs01_.addZOT( valueOfFunc(Psi, [](auto psi) { return -3.0 * Math::pow<2>(psi); }, 2) );
    opLhs01_.addSOT( 2.0 );
    auto opLhs01 = std::make_pair(opLhs01_, true);

    // < -2*psi^3, v >
    auto opRhs0 = Op::zot(valueOfFunc(Psi, [](auto psi) { return -2.0 * Math::pow<3>(psi); }, 3));

    // < psi, v >
    auto opRhs1 = Op::zot(valueOf(Psi));

    // < u, v >
    auto opM = std::make_pair(Op::zot(1.0), false);

    // < grad(u), grad(v) >
    auto opL = std::make_pair(Op::sot(1.0), false);

    // < M0*tau * grad(u), grad(v) >
    double* tauPtr = adaptInfo.getTimestepPtr();
    auto opL0 = std::make_pair(Op::sot([M0, tauPtr](auto /*x*/) { return M0 * (*tauPtr); }), false);


    // === Define the linear system ===

    prob.addMatrixOperator({{ {{0,0}, opM},  {{0,1}, opLhs01}, {{0,2}, opL},
                              {{1,0}, opL0}, {{1,1}, opM},
                                             {{2,1}, opL},     {{2,2}, opM}}});

    prob.addVectorOperator({{ {0, opRhs0},
                              {1, opRhs1} }});

    // === set the initial solution ===
    Mu = 0.0;
    Nu = 0.0;
    std::srand( 12345 );
    Psi.interpol([psi_mean](auto const& x) { return 0.5*(std::rand()/double(RAND_MAX) - 0.5) + psi_mean; });

    // === configure preconditioner, if necessary ===
    if (creator->precon) {
      creator->precon->setData(tauPtr, M0);
    }

    AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo);
    adapt.adapt();

    msg("|mu|  = ", two_norm(Mu.getVector()), ", ",
        "|nu|  = ", two_norm(Nu.getVector()), ", ",
        "|psi| = ", two_norm(Psi.getVector()));

    AMDiS::finalize();
    return 0;
}