Skip to content
Snippets Groups Projects
Commit a58f1ace authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Remove laplace tests

We are in a state, where Poisson is totally sufficient as a minimal problem.
parent 21735855
No related branches found
No related tags found
No related merge requests found
dune_add_formcompiler_system_test(UFLFILE laplace.ufl
BASENAME laplace
INIFILE laplace.mini)
dune_add_formcompiler_system_test(UFLFILE laplace_dg.ufl
BASENAME laplace_dg
INIFILE laplace_dg.mini)
add_executable(laplace_dg_ref reference_main.cc)
set_target_properties(laplace_dg_ref PROPERTIES EXCLUDE_FROM_ALL 1)
__name = laplace_{__exec_suffix}
__exec_suffix = symdiff, numdiff | expand num
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 4 4
elementType = simplical
printmatrix = true
[formcompiler.r]
numerical_jacobian = 0, 1 | expand num
V = FiniteElement("CG", "triangle", 1)
u = TrialFunction(V)
v = TestFunction(V)
r = inner(grad(u), grad(v))*dx
__name = laplace_dg_{__exec_suffix}
__exec_suffix = numdiff, symdiff | expand num
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 2 2
elementType = simplical
printmatrix = true
[formcompiler.r]
numerical_jacobian = 1, 0 | expand num
cell = triangle
V = FiniteElement("DG", cell, 1)
u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(cell)('+')
# penalty factor
gamma = 1.0
# SIPG: -1.0, IIPG: 0.0, NIPG: 1.0
theta = 1.0
r = inner(grad(u), grad(v))*dx \
+ inner(n, avg(grad(u)))*jump(v)*dS \
+ gamma*jump(u)*jump(v)*dS \
- theta*jump(u)*inner(avg(grad(v)), n)*dS \
- inner(n, grad(u))*v*ds \
+ gamma*u*v*ds \
+ theta*u*inner(grad(v), n)*ds
forms = [r]
#ifndef _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
#define _HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
#include <dune/pdelab/gridoperator/gridoperator.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/backend/istl.hh>
#include <dune/pdelab/gridfunctionspace/vtk.hh>
#include <dune/grid/uggrid.hh>
#include <dune/pdelab/backend/istl.hh>
#include <string>
#include <dune/common/parametertree.hh>
#include <dune/pdelab/finiteelementmap/opbfem.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/pdelab/stationary/linearproblem.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/testtools/gridconstruction.hh>
#include <dune/pdelab/localoperator/convectiondiffusiondg.hh>
#include <dune/pdelab/localoperator/convectiondiffusionparameter.hh>
template<typename GV, typename RF>
class CDProb
{
typedef Dune::PDELab::ConvectionDiffusionBoundaryConditions::Type BCType;
public:
typedef Dune::PDELab::ConvectionDiffusionParameterTraits<GV,RF> Traits;
//! tensor diffusion coefficient
typename Traits::PermTensorType
A (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
{
typename Traits::PermTensorType I;
for (std::size_t i=0; i<Traits::dimDomain; i++)
for (std::size_t j=0; j<Traits::dimDomain; j++)
I[i][j] = (i==j) ? 1 : 0;
return I;
}
//! velocity field
typename Traits::RangeType
b (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
{
typename Traits::RangeType v(0.0);
return v;
}
//! sink term
typename Traits::RangeFieldType
c (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
{
return 0.0;
}
//! source term
typename Traits::RangeFieldType
f (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
{
return 0.0;
}
//! boundary condition type function
BCType
bctype (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
{
return Dune::PDELab::ConvectionDiffusionBoundaryConditions::Dirichlet;
}
//! Dirichlet boundary condition value
typename Traits::RangeFieldType
g (const typename Traits::ElementType& e, const typename Traits::DomainType& x) const
{
return 0.0;
}
//! Neumann boundary condition
typename Traits::RangeFieldType
j (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
{
return 0.0;
}
//! outflow boundary condition
typename Traits::RangeFieldType
o (const typename Traits::IntersectionType& is, const typename Traits::IntersectionDomainType& x) const
{
return 0.0;
}
};
void driver(int argc, char** argv){ typedef Dune::PDELab::ISTL::VectorBackend<Dune::PDELab::ISTL::Blocking::none, 1> VectorBackend;
static const int dim = 2;
typedef Dune::UGGrid<dim> Grid;
typedef Grid::LeafGridView GV;
typedef Grid::ctype DF;
typedef double R;
typedef Dune::PDELab::OPBLocalFiniteElementMap<DF, R, 1, dim, Dune::GeometryType::simplex> DG1_FEM;
typedef Dune::PDELab::NoConstraints NoConstraintsAssembler;
typedef Dune::PDELab::GridFunctionSpace<GV, DG1_FEM, NoConstraintsAssembler, VectorBackend> DG1_DIRICHLET_GFS;
Dune::ParameterTree initree;
Dune::ParameterTreeParser::readINITree(argv[1], initree);
IniGridFactory<Grid> factory(initree);
std::shared_ptr<Grid> grid = factory.getGrid();
GV gv = grid->leafGridView();
DG1_FEM dg1_fem;
DG1_DIRICHLET_GFS dg1_dirichlet_gfs(gv, dg1_fem);
dg1_dirichlet_gfs.name("bla");
typedef Dune::SubsamplingVTKWriter<GV> VTKWriter;
int sublevel = initree.get<int>("vtk.subsamplinglevel", 0);
VTKWriter vtkwriter(gv, sublevel);
using LOP = Dune::PDELab::ConvectionDiffusionDG<CDProb<GV, R>, DG1_FEM>;
typedef DG1_DIRICHLET_GFS::ConstraintsContainer<R>::Type DG1_CC;
typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MatrixBackend;
typedef Dune::PDELab::GridOperator<DG1_DIRICHLET_GFS, DG1_DIRICHLET_GFS, LOP, MatrixBackend, DF, R, R, DG1_CC, DG1_CC> GO;
typedef GO::Traits::Domain V;
V x(dg1_dirichlet_gfs);
x = 0.0;
std::string vtkfile = initree.get<std::string>("wrapper.vtkcompare.name", "output");
typedef Dune::PDELab::ISTLBackend_SEQ_UMFPack LinearSolver;
typedef Dune::PDELab::StationaryLinearProblemSolver<GO, LinearSolver, V> SLP;
DG1_CC dg1_cc;
dg1_cc.clear();
CDProb<GV, R> params;
LOP lop(params, Dune::PDELab::ConvectionDiffusionDGMethod::SIPG);
int generic_dof_estimate = 6 * dg1_dirichlet_gfs.maxLocalSize();
int dofestimate = initree.get<int>("istl.number_of_nnz", generic_dof_estimate);
MatrixBackend mb(dofestimate);
GO go(dg1_dirichlet_gfs, dg1_cc, dg1_dirichlet_gfs, dg1_cc, lop, mb);
std::cout << "gfs with " << dg1_dirichlet_gfs.size() << " dofs generated "<< std::endl;
std::cout << "cc with " << dg1_cc.size() << " dofs generated "<< std::endl;
LinearSolver ls(false);
double reduction = initree.get<double>("reduction", 1e-12);
SLP slp(go, ls, x, reduction);
slp.apply();
typedef typename GO::Traits::Jacobian M;
M m(go);
go.jacobian(x,m);
using Dune::PDELab::Backend::native;
Dune::printmatrix(std::cout, native(m),"global stiffness matrix","row",9,1);
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, dg1_dirichlet_gfs, x);
vtkwriter.write(vtkfile, Dune::VTK::ascii);
}
#endif //_HOME_DOMINIC_DUNE_BUILD_DUNE_PERFTOOL_TEST_LAPLACE_LAPLACE_DG_SYMDIFF_DRIVER_HH
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/exceptions.hh>
#include"/home/dominic/dune/dune-perftool/test/laplace/reference_driver.hh"
int main(int argc, char** argv)
{
try{
//Maybe initialize Mpi
Dune::MPIHelper& helper = Dune::MPIHelper::instance(argc, argv);
if(Dune::MPIHelper::isFake)
std::cout<< "This is a sequential program." << std::endl;
else
std::cout<<"I am rank "<<helper.rank()<<" of "<<helper.size()
<<" processes!"<<std::endl;
driver(argc, argv);
return 0;
}
catch (Dune::Exception &e){
std::cerr << "Dune reported error: " << e << std::endl;
return 1;
}
catch (...){
std::cerr << "Unknown exception thrown!" << std::endl;
return 1;
}
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment