Commit 4594dccb authored by Andreas Nüßing's avatar Andreas Nüßing
Browse files

add finite volume to udg prolongation

parent 10cffbcf
......@@ -8,6 +8,7 @@
#include <duneuro/common/dg_solver_backend_factory.hh>
#include <duneuro/common/edge_norm_provider.hh>
#include <duneuro/common/flags.hh>
#include <duneuro/common/fv_space.hh>
#include <duneuro/common/random.hh>
#include <duneuro/common/thread_safe_linear_problem_solver.hh>
#include <duneuro/io/data_tree.hh>
......@@ -45,6 +46,9 @@ namespace duneuro
using EdgeNormProvider = MultiEdgeNormProvider;
using LocalOperator = ConvectionDiffusion_DG_LocalOperator<Problem, EdgeNormProvider>;
using Assembler = GalerkinGlobalAssembler<FunctionSpace, LocalOperator, DF, RF, JF>;
using FVSpace =
duneuro::FVSpace<typename VC::GridView, BasicTypeFromElementType<elementType>::value,
Dune::SolverCategory::sequential>;
using FirstOrderSpace =
CGFirstOrderSpace<typename VC::GridView, BasicTypeFromElementType<elementType>::value,
Dune::SolverCategory::sequential>;
......@@ -72,7 +76,11 @@ namespace duneuro
false, config.get<unsigned int>("intorderadd", 0))
, assembler_(functionSpace_, localOperator_, 2 * VC::dim + 1)
, firstOrderSpace_(volumeConductor_->grid(), volumeConductor_->gridView())
, solverBackend_(make_dg_solver_backend(*assembler_, firstOrderSpace_.getGFS(), config))
, fvSpace_(volumeConductor_->grid(), volumeConductor_->gridView())
, solverBackend_(
config.get<std::string>("prolongation_type", "linear") == "linear" ?
make_dg_solver_backend(*assembler_, firstOrderSpace_.getGFS(), config) :
make_dg_solver_backend(*assembler_, fvSpace_.getGFS(), config))
, linearSolver_(*assembler_, config)
{
assert(volumeConductor_);
......@@ -105,6 +113,7 @@ namespace duneuro
typename Traits::LocalOperator localOperator_;
typename Traits::Assembler assembler_;
typename Traits::FirstOrderSpace firstOrderSpace_;
typename Traits::FVSpace fvSpace_;
std::unique_ptr<SolverBackendInterface<typename Traits::Assembler::GO::Traits::Jacobian,
typename Traits::DomainDOFVector,
typename Traits::RangeDOFVector>>
......
......@@ -3,6 +3,7 @@
#include <dune/common/std/memory.hh>
#include <duneuro/common/fv_to_dg_prolongation.hh>
#include <duneuro/common/solver_backends.hh>
namespace duneuro
......@@ -19,14 +20,30 @@ namespace duneuro
if (solverType == "cg") {
auto precoditionerType = config.get<std::string>("preconditioner_type", "amg");
if (precoditionerType == "amg") {
auto dgSmootherType = config.get<std::string>("dg_smoother_type", "ssor");
if (dgSmootherType == "ssor") {
return Dune::Std::make_unique<AMG4DGSolverBackend<DGGO, CGGFS,
Dune::PDELab::CG2DGProlongation,
Dune::SeqSSOR, Dune::CGSolver>>(
dggo, cggfs, config);
auto prolongationType = config.get<std::string>("prolongation_type", "linear");
if (prolongationType == "linear") {
auto dgSmootherType = config.get<std::string>("dg_smoother_type", "ssor");
if (dgSmootherType == "ssor") {
return Dune::Std::make_unique<AMG4DGSolverBackend<DGGO, CGGFS,
Dune::PDELab::CG2DGProlongation,
Dune::SeqSSOR, Dune::CGSolver>>(
dggo, cggfs, config);
} else {
DUNE_THROW(Dune::Exception, "unknown dg smoother type \"" << dgSmootherType << "\"");
}
} else if (prolongationType == "fv") {
auto dgSmootherType = config.get<std::string>("dg_smoother_type", "ssor");
if (dgSmootherType == "ssor") {
return Dune::Std::make_unique<AMG4DGSolverBackend<DGGO, CGGFS,
duneuro::FV2DGProlongation,
Dune::SeqSSOR, Dune::CGSolver>>(
dggo, cggfs, config);
} else {
DUNE_THROW(Dune::Exception, "unknown dg smoother type \"" << dgSmootherType << "\"");
}
} else {
DUNE_THROW(Dune::Exception, "unknown dg smoother type \"" << dgSmootherType << "\"");
DUNE_THROW(Dune::Exception, "unknown prolongation type type \"" << prolongationType
<< "\"");
}
} else if (precoditionerType == "ssor") {
return Dune::Std::make_unique<CGSSORSolverBackend<M, D, R>>(config);
......
......@@ -16,6 +16,7 @@
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <duneuro/common/cg_to_udg_prolongation.hh>
#include <duneuro/common/fv_to_udg_prolongation.hh>
namespace duneuro
{
......@@ -48,6 +49,20 @@ namespace duneuro
Dune::printmatrix(out, matrix, "jacobian", "r");
}
}
enum class UDGProlongationType { linear, fv };
static inline UDGProlongationType udg_prolongation_from_string(const std::string& value)
{
if (value == "linear") {
return UDGProlongationType::linear;
} else if (value == "fv") {
return UDGProlongationType::fv;
} else {
DUNE_THROW(Dune::Exception, "unknown prolongation type \"" << value << "\"");
}
}
/** Sequential solver backend for using AMG for DG in PDELab
The template parameters are:
......@@ -101,6 +116,8 @@ namespace duneuro
unsigned int solverIterations;
bool additive;
UDGProlongationType prolongationType;
bool reuse;
bool firstapply;
......@@ -122,6 +139,8 @@ namespace duneuro
, dgSmootherRelaxation(params.get<field_type>("dg_smoother_relaxation", 1.0))
, solverIterations(params.get<unsigned int>("max_iterations", 5000))
, additive(params.get<bool>("additive", false))
, prolongationType(
udg_prolongation_from_string(params.get<std::string>("prolongation_type")))
, reuse(true)
, firstapply(true)
{
......@@ -144,7 +163,7 @@ namespace duneuro
if (config.hasKey("additive"))
additive = config.get<bool>("additive");
if (config.hasKey("verbose"))
additive = config.get<int>("verbose");
verbose = config.get<int>("verbose");
}
/*! \brief compute global norm of a vector
......@@ -186,10 +205,14 @@ namespace duneuro
// no need to set acg here back to zero, this is done in matMultmat
if (reuse == false || firstapply == true) {
// pmatrix = constructCG2UDGProlongation(native(A));
pmatrix = compute_cg_to_udg_prolongation(udggfs, st);
//auto pmatrix_new = compute_cg_to_udg_prolongation_new(udggfs, st);
if (prolongationType == UDGProlongationType::linear) {
pmatrix = compute_cg_to_udg_prolongation(udggfs, st);
} else {
pmatrix = compute_fv_to_udg_prolongation(udggfs, st);
}
// auto pmatrix_new = compute_cg_to_udg_prolongation_new(udggfs, st);
//*pmatrix_new -= *pmatrix;
//std::cout << "difference: " << pmatrix_new->frobenius_norm() << "\n";
// std::cout << "difference: " << pmatrix_new->frobenius_norm() << "\n";
// Dune::writeMatrixToMatlab(*pmatrix, "pmatrix.mat");
// Dune::writeMatrixToMatlab(native(A), "A.mat");
std::cout << "A size: " << native(A).N() << "x" << native(A).M() << std::endl;
......
#ifndef DUNEURO_UDG_MULTI_PHASE_SPACE_HH
#define DUNEURO_UDG_MULTI_PHASE_SPACE_HH
#include <dune/localfunctions/lagrange/p0.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/finiteelementmap/l2orthonormal.hh>
#include <dune/pdelab/finiteelementmap/qkdg.hh>
......@@ -67,6 +69,7 @@ namespace duneuro
UDGQkMultiPhaseSpace(const UDGQkMultiPhaseSpace&) = delete;
UDGQkMultiPhaseSpace& operator=(const UDGQkMultiPhaseSpace&) = delete;
private:
GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_;
......@@ -76,5 +79,68 @@ namespace duneuro
std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_;
std::shared_ptr<GFS> gfs_;
};
template <class TGV, class N, int phases>
class UDGP0MultiPhaseSpace
{
public:
typedef TGV GV;
enum { dim = GV::dimension };
typedef typename GV::ctype ctype;
typedef N NT;
typedef Dune::P0LocalFiniteElement<ctype, NT, dim> LFE;
typedef VirtualSubTriangulation<GV> SubTriangulation;
enum { blockSize = 1 };
typedef Dune::PDELab::istl::VectorBackend<Dune::PDELab::istl::Blocking::fixed, blockSize> VBE;
typedef Dune::PDELab::UnfittedFiniteElementMapTraits<LFE, typename SubTriangulation::EntityPart>
UFEMTraits;
typedef Dune::PDELab::UnfittedFiniteElementMap<UFEMTraits, SubTriangulation> FEM;
typedef Dune::PDELab::LeafOrderingTag<UDGLeafOrderingParams> LeafOrderingTag;
typedef Dune::PDELab::GridFunctionSpace<GV, FEM, Dune::PDELab::NoConstraints, VBE,
LeafOrderingTag>
DomainGFS;
typedef Dune::PDELab::istl::VectorBackend<> PVBE;
typedef Dune::PDELab::ordering::Chunked<Dune::PDELab::EntityBlockedOrderingTag> OrderingTag;
typedef Dune::PDELab::PowerGridFunctionSpace<DomainGFS, phases, PVBE, OrderingTag> GFS;
typedef typename Dune::PDELab::BackendVectorSelector<GFS, N>::Type DOF;
UDGP0MultiPhaseSpace(const GV& gv, std::shared_ptr<const SubTriangulation> subTriangulation)
: gridView_(gv)
, entitySet_(gridView_)
, subTriangulation_(subTriangulation)
, lfe_(gv.template begin<0>()->geometry().type())
{
for (unsigned int i = 0; i < phases; ++i) {
fems_[i] = std::make_shared<FEM>(lfe_, *subTriangulation_, i);
domainGfss_[i] = std::make_shared<DomainGFS>(entitySet_, *(fems_[i]));
}
gfs_ = make_power_gfs(domainGfss_, PVBE(), OrderingTag(blockSize));
gfs_->ordering();
}
// return gfs reference
GFS& getGFS()
{
return *gfs_;
}
// return gfs reference const version
const GFS& getGFS() const
{
return *gfs_;
}
UDGP0MultiPhaseSpace(const UDGP0MultiPhaseSpace&) = delete;
UDGP0MultiPhaseSpace& operator=(const UDGP0MultiPhaseSpace&) = delete;
private:
GV gridView_;
Dune::PDELab::AllEntitySet<GV> entitySet_;
std::shared_ptr<const SubTriangulation> subTriangulation_;
LFE lfe_;
std::array<std::shared_ptr<FEM>, phases> fems_;
std::array<std::shared_ptr<DomainGFS>, phases> domainGfss_;
std::shared_ptr<GFS> gfs_;
};
}
#endif // DUNEURO_UDG_MULTI_PHASE_SPACE_HH
Markdown is supported
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