Skip to content
GitLab
Commit 048ac58b authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Merge branch 'issue/remove_old_examples' into 'develop' 

old examples that do not compile any more removed

See merge request spraetor/dune-amdis!42
parents 7c4aa5c5 c61b2e0d
Hide whitespace changes
Inline Side-by-side
examples/old/BlockPreconditioner.hpp deleted 100644 → 0
View file @ 7c4aa5c5
#pragma once
#include <amdis/AMDiS.hpp>
#include <amdis/LinearAlgebra.hpp>
namespace AMDiS
{
/// Basis preconditioner structure for block-preconditioners
template <class Matrix, class Vector>
class BlockPreconditioner
: public PreconditionerInterface<Matrix, Vector>
{
public:
using Self = BlockPreconditioner;
using Super = PreconditionerInterface<Matrix, Vector>;
using PreconBase = Super;
public:
BlockPreconditioner() = default;
/// extract iranges from BlockMatrix to be used to extract sub-vectors and sub-matrices.
virtual void init(Matrix const& A_) override
{
A = &A_;
A_.getRanges(rows, cols);
initialized = true;
}
virtual void solve(Vector const& b, Vector& x) const override
{
error_exit("Must be implemented in derived classes!\n");
}
virtual void adjoint_solve(Vector const& x, Vector& y) const override
{
error_exit("Must be implemented in derived classes!\n");
}
template <size_t I>
mtl::irange const& getRowRange(const index_t<I> i = {}) const
{
test_exit_dbg(initialized,
"BlockPreconditioner not initialized. Call init() before.");
return std::get<I>(rows);
}
template <size_t I>
mtl::irange const& getColRange(const index_t<I> i = {}) const
{
test_exit_dbg(initialized,
"BlockPreconditioner not initialized. Call init() before.");
return std::get<I>(cols);
}
protected:
Matrix const* A = NULL;
bool initialized = false;
std::array<mtl::irange, Matrix::N()> rows;
std::array<mtl::irange, Matrix::M()> cols;
};
} // end namespace AMDiS
examples/old/MTLPfcPrecon.hpp deleted 100644 → 0
View file @ 7c4aa5c5
#pragma once
#include "BlockPreconditioner.hpp"
#include <amdis/common/Literals.hpp>
namespace AMDiS
{
template <class Matrix, class Vector>
class MTLPfcPrecon
: public BlockPreconditioner<Matrix, Vector>
{
public:
using Self = MTLPfcPrecon;
using Super = BlockPreconditioner<Matrix, Vector>;
using PreconBase = typename Super::PreconBase;
using MTLMatrix = BaseMatrix_t<Matrix>;
using MTLVector = Vector;
using LinearSolverType = LinearSolverInterface<MTLMatrix, MTLVector>;
public:
struct Creator : CreatorInterface<PreconBase>
{
virtual std::shared_ptr<PreconBase> create() override
{
precon = std::make_shared<Self>();
return precon;
}
std::shared_ptr<Self> precon = NULL;
};
public:
MTLPfcPrecon()
: Super()
{
std::string solverNameM = "cg",
solverNameMpL = "cg",
solverNameMpL2 = "cg";
Parameters::get("precon_pfc->M->solver", solverNameM);
Parameters::get("precon_pfc->MpL->solver", solverNameMpL);
Parameters::get("precon_pfc->MpL2->solver", solverNameMpL2);
CreatorInterfaceName<LinearSolverType>* creator;
creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameM, "precon_pfc->M->solver") );
solverM = creator->create("precon_pfc->M");
creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL, "precon_pfc->MpL->solver") );
solverMpL = creator->create("precon_pfc->MpL");
creator = named( CreatorMap<LinearSolverType>::getCreator(solverNameMpL2, "precon_pfc->MpL2->solver") );
solverMpL2 = creator->create("precon_pfc->MpL2");
}
virtual ~MTLPfcPrecon()
{
exit();
}
void setData(double* tau_, double M0_ = 1.0)
{
tau = tau_;
M0 = M0_;
}
// implementation in MTLPfcPrecon.inc.hpp
virtual void init(Matrix const& A) override;
virtual void exit() override {}
// implementation in MTLPfcPrecon.inc.hpp
virtual void solve(MTLVector const& b, MTLVector& x) const override;
template <typename VectorX, typename VectorB>
void solveM(VectorX& x, const VectorB& b) const
{
SolverInfo solverInfo("precon_pfc->M");
solverM->solve(getM(), x, b, solverInfo);
}
template <typename VectorX, typename VectorB>
void solveMpL(VectorX& x, const VectorB& b) const
{
SolverInfo solverInfo("precon_pfc->MpL");
solverMpL->solve(MpL, x, b, solverInfo);
}
template <typename VectorX, typename VectorB>
void solveMpL2(VectorX& x, const VectorB& b) const
{
SolverInfo solverInfo("precon_pfc->MpL2");
solverMpL2->solve(MpL2, x, b, solverInfo);
}
MTLMatrix const& getM() const { return matM ? *matM : (*Super::A)(_2, _2); } // < u, v >
MTLMatrix const& getL0() const { return matL0 ? *matL0 : (*Super::A)(_1, _0); } // < M0*tau*grad(u), grad(v) >
MTLMatrix const& getL() const { return matL ? *matL : (*Super::A)(_2, _1); } // < grad(u), grad(v) >
double getTau() const { return *tau; }
protected:
double* tau = NULL;
double M0 = 1.0;
MTLMatrix* matM = NULL;
MTLMatrix* matL = NULL;
MTLMatrix* matL0 = NULL;
MTLMatrix MpL;
MTLMatrix MpL2;
mutable MTLVector y0;
mutable MTLVector y1;
mutable MTLVector tmp;
std::shared_ptr<LinearSolverType> solverM;
std::shared_ptr<LinearSolverType> solverMpL;
std::shared_ptr<LinearSolverType> solverMpL2;
};
} // end namespace AMDiS
#include "MTLPfcPrecon.inc.hpp"
examples/old/MTLPfcPrecon.inc.hpp deleted 100644 → 0
View file @ 7c4aa5c5
namespace AMDiS
{
template <class Matrix, class Vector>
void MTLPfcPrecon<Matrix, Vector>::init(Matrix const& A)
{
test_exit(tau, "Preconditioner not initialized!");
Super::init(A);
using size_type = typename MTLMatrix::size_type;
size_type n = num_rows(getM());
double delta = std::sqrt(M0 * getTau());
msg("delta = {}", delta);
// helper-matrix MpL = M + delta*L
MpL.change_dim(n, n);
MpL = getM() + (1.0/delta) * getL0();
// helper-matrix MpL = M + sqrt(delta)*L
MpL2.change_dim(n, n);
MpL2 = getM() + std::sqrt(delta) * getL();
// temporary variables
y0.change_dim(num_rows(getM()));
y1.change_dim(num_rows(getM()));
tmp.change_dim(num_rows(getM()));
}
template <class Matrix, class Vector>
void MTLPfcPrecon<Matrix, Vector>::solve(Vector const& b, Vector& x) const
{
x.change_dim(num_rows(b));
Vector x0(x[Super::getColRange(0_c)]);
Vector x1(x[Super::getColRange(1_c)]);
Vector x2(x[Super::getColRange(2_c)]);
const Vector b0(b[Super::getRowRange(0_c)]);
const Vector b1(b[Super::getRowRange(1_c)]);
const Vector b2(b[Super::getRowRange(2_c)]);
double delta = std::sqrt(M0 * getTau());
solveM(y0, b0); // M*y0 = b0
y1 = getL0() * y0; // y1 = K*y0
tmp = b1 - y1; // tmp := b1 - tau*y1
solveMpL(y1, tmp); // (M + delta*K) * y1 = tmp
x0 = y0 + (1.0/delta)*y1; // x0 = y0 + (1/delta)*y1
tmp = getM() * y1; // tmp := M*y1
solveMpL2(y1, tmp); // (M+eps*sqrt(delta)K) * y1 = tmp
tmp = getM() * y1; // tmp := M*y1
solveMpL2(x1, tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
x0-= (1.0/delta)*x1; // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
tmp = b2 - getL() * x1; // tmp := b2 - K*x1
solveM(x2, tmp);
}
} // end namespace AMDiS
examples/old/PfcPrecon.hpp deleted 100644 → 0
View file @ 7c4aa5c5
#pragma once
#include <tuple>
#include <cmath>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/multitypeblockmatrix.hh>
#include <dune/istl/multitypeblockvector.hh>
#include <dune/istl/preconditioner.hh>
#include <amdis/common/Mpl.hpp>
namespace AMDiS {
template <class Matrix, class Vector>
class PfcPrecon
: public Dune::Preconditioner<Vector, Vector>
{
using SubMatrix = std::decay_t<decltype( std::declval<Matrix>()[_0][_0] )>;
using SubVector = std::decay_t<decltype( std::declval<Vector>()[_0] )>;
public:
PfcPrecon(Matrix const& matrix, double* tauPtr, double M0)
: matrix(matrix)
, tauPtr(tauPtr)
, M0(M0)
, matL0(matrix[_1][_0])
, matL (matrix[_2][_1])
, matM (matrix[_2][_2])
{}
virtual void pre(Vector& x, Vector& b) override
{
double delta = std::sqrt(M0 * (*tauPtr));
matMpL = matM;
matMpL.axpy(1.0/delta, matL0); // => MpL = M + 1/delta * L0
matMpL2 = matM;
matMpL2.axpy(std::sqrt(delta), matL);
y0.resize(matM.N());
y1.resize(matM.N());
tmp.resize(matM.N());
}
virtual void apply(Vector& x, const Vector& b) override
{
double delta = std::sqrt(M0 * (*tauPtr));
solve(matM, y0, b[_0]); // M*y0 = b0
matL0.mv(y0, y1); // y1 = K*y0
tmp = b[_1];
tmp-= y1; // tmp := b1 - tau*y1
solve(matMpL, y1, tmp); // (M + delta*K) * y1 = tmp
x[_0] = y0;
x[_0].axpy(1.0/delta, y1); // x0 = y0 + (1/delta)*y1
matM.mv(y1, tmp); // tmp := M*y1
solve(matMpL2, y1, tmp); // (M+eps*sqrt(delta)K) * y1 = tmp
matM.mv(y1, tmp); // tmp := M*y1
solve(matMpL2, x[_1], tmp); // (M+eps*sqrt(delta)K) * x1 = tmp
x[_0].axpy(-1.0/delta, x[_1]); // x0 = x0 - (1/delta)*x1 = y0 + (1/delta)*(y1 - x1)
matL.mv(x[_1], y1);
tmp = b[_2];
tmp-= y1; // tmp := b2 - K*x1
solve(matM, x[_2], tmp);
}
virtual void post(Vector& x) override
{
}
private:
template <class Mat, class Vec>
void solve(Mat const& A, Vec& x, Vec b)
{
Dune::MatrixAdapter<Mat, Vec, Vec> op(A);
Dune::SeqJac<Mat, Vec, Vec> precon(A, 1, 1.0);
Dune::CGSolver<Vec> solver(op, precon, 1.e-3, 20, 0);
Dune::InverseOperatorResult statistics;
solver.apply(x, b, statistics);
}
private:
Matrix const& matrix;
double* tauPtr;
double M0;
SubMatrix const& matL0;
SubMatrix const& matL;
SubMatrix const& matM;
SubMatrix matMpL;
SubMatrix matMpL2;
SubVector y0, y1, tmp;
};
} // end namespace AMDiS
examples/old/pfc.cc deleted 100644 → 0
View file @ 7c4aa5c5
// -*- 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| = {}, |nu| = {}, |psi| = {}",
two_norm(Mu.getVector()), two_norm(Nu.getVector()), two_norm(Psi.getVector()));
AMDiS::finalize();
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment