Commit 1114932d authored by Peter Bastian's avatar Peter Bastian
Browse files

the geneo example

parent d4805639
......@@ -3,4 +3,4 @@ target_link_dune_default_libraries(example)
add_executable(example_geneo example_geneo.cc)
target_link_dune_default_libraries(example_geneo)
dune_symlink_to_source_files(FILES example.yml)
dune_symlink_to_source_files(FILES example.yml example.ini)
[geneo]
coordinate_partitioning=true
parts=4 1
subdomains= 4 1
overlap=2
coarse_overlap=0
pum=distance
fineGEVPrhs=pu
coarseGEVPrhs=pu
extensionmethod=vertex
view=0
vizbasisvectors=false
arpack_tolerance = 1e-8
n_eigenvectors_fine_computed = 40
n_eigenvectors_fine_used = 25
n_eigenvectors_coarse_computed = 120 # eigen computes all anyway
n_eigenvectors_coarse_used = 24
verbose = 0
abs_zero_ker = 0.0
regularization_ker = 1e-8
merge_disconnected=false
drop_small_overlap=1
coarseeigensolver=straightarpack
eigenvalue_fine_threshold=0.15
eigenvalue_coarse_threshold=0.3
cycle=additive
[solver]
itmax=1000
reduction=1e-13
type=geneo
......@@ -6,8 +6,8 @@ grid:
- 4.0
- 1.0
N:
- 40
- 10
- 80
- 20
refinement: 0
material:
-
......@@ -35,11 +35,6 @@ reinforced_operator:
- -0.5
- 0.749
youngs_modulus: 10000
stabilization_parameter: 1.0
geneo:
subdomains:
- 4
- 1
overlap: 2
extensionmethod: vertex
drop_small_overlap: 1
stabilization_parameter: 100.0
enable_skeleton: true
\ No newline at end of file
......@@ -12,6 +12,7 @@
#include<dune/fiber-elasticity/vonmises.hh>
// General Dune stuff
#include<dune/common/parametertreeparser.hh>
#include<dune/common/parallel/mpihelper.hh>
#include<dune/grid/io/file/gmshreader.hh>
#include<dune/grid/uggrid.hh>
......@@ -100,25 +101,28 @@ int main(int argc, char** argv)
lop.setCoefficientTraction(Dune::stackobject_to_shared_ptr(vgfs), Dune::stackobject_to_shared_ptr(traction_force));
lop.setCoefficientForce(Dune::stackobject_to_shared_ptr(vgfs), Dune::stackobject_to_shared_ptr(body_force));
// Linear Solver Setup
using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_UMFPack;
using StationaryLinearProblemSolver = Dune::PDELab::StationaryLinearProblemSolver<GridOperator, LinearSolver, Vector>;
LinearSolver ls(0);
StationaryLinearProblemSolver slp(go, ls, displacement, 1e-12);
//==================
// Begin Geneo stuff
//==================
// read parameters
std::vector<size_t> subdomains = config["geneo"]["subdomains"].template as<std::vector<size_t>>(std::vector<size_t>(2,1));
std::cout << "subdomains: "; for (auto p : subdomains) std::cout << p << " "; std::cout << std::endl;
size_t overlap = config["geneo"]["overlap"].template as<size_t>(1);
std::string extensionmethod = config["geneo"]["extensionmethod"].template as<std::string>(std::string("vertex"));
size_t drop_small_overlap = config["geneo"]["drop_small_overlap"].template as<size_t>(1);
Dune::ParameterTree ptree;
Dune::ParameterTreeParser ptreeparser;
ptreeparser.readINITree("example.ini",ptree);
ptreeparser.readOptions(argc,argv,ptree);
std::vector<size_t> subdomains = ptree.get<std::vector<size_t>>("geneo.subdomains");
size_t overlap = ptree.get("geneo.overlap",(size_t)1);
std::string extensionmethod = ptree.get("geneo.extensionmethod","vertex");
size_t drop_small_overlap = ptree.get("geneo.drop_small_overlap",(size_t)1);
bool coordinate_partitioning = ptree.get("geneo.coordinate_partitioning",(bool)false);
std::vector<size_t> parts = ptree.get<std::vector<size_t>>("geneo.parts");
// load balancing and subdomain information for geneo
auto pdd = std::make_shared<DomainDecomposition>(gv,Dune::MPIHelper::getLocalCommunicator(),subdomains[0],overlap,extensionmethod,drop_small_overlap,false);
std::shared_ptr<DomainDecomposition> pdd;
if (coordinate_partitioning)
pdd = std::make_shared<DomainDecomposition>(gv,subdomains[0],parts,overlap,extensionmethod,drop_small_overlap,false);
else
pdd = std::make_shared<DomainDecomposition>(gv,Dune::MPIHelper::getLocalCommunicator(),subdomains[0],overlap,extensionmethod,drop_small_overlap,false);
auto pdddm = std::make_shared<DomainDecompositionDOFMapper>(vgfs,pdd);
// set up linear system in PDELab components
......@@ -158,6 +162,7 @@ int main(int argc, char** argv)
wgop.jacobian(displacement,A); // assembles the global stiffness matrix and the snippets
auto assemble_time = timer.elapsed();
std::cout << "SNIPPET ASSEMBLE TIME=" << assemble_time << std::endl;
//Dune::printmatrix(std::cout, istlA, " A", "");
// global vector with dirichlet boundary conditions
// contains 0.0 for dof on Dirichlet boundary, otherwise 1.0
......@@ -180,35 +185,102 @@ int main(int argc, char** argv)
if ((*pistldirichlet_mask)[cIt.index()][compj]==0.0)
(*cIt)[compi][compj] = 0.0;
}
bool sym = is_symmetric(istlA,"A",true,1e-12);
bool sym = is_symmetric(istlA,"A",false,1e-8);
std::cout << "matrix is symmetric: " << sym << std::endl;
std::cout << "matrix norm: " << istlA.infinity_norm() << std::endl;
// for (size_t i=0; i<istlA.N(); i++)
// {
// auto cIt = istlA[i].begin();
// auto cEndIt = istlA[i].end();
// for (; cIt!=cEndIt; ++cIt)
// if (cIt.index()==i)
// for (size_t compi=0; compi<blocksize; compi++)
// std::cout << i << " " << compi << " " << (*cIt)[compi][compi] << std::endl;
// }
size_t coarse_overlap = ptree.get("geneo.coarse_overlap",(size_t)0);
std::string pum = ptree.get("geneo.pum","standard");
std::string fineGEVPrhs = ptree.get("geneo.fineGEVPrhs","pu");
std::string coarseGEVPrhs = ptree.get("geneo.coarseGEVPrhs","pu");
std::string coarseeigensolver = ptree.get("geneo.coarseeigensolver","eigen");
std::string cycle = ptree.get("geneo.cycle","additive");
size_t n_eigenvectors_fine_computed = ptree.get("geneo.n_eigenvectors_fine_computed",(size_t)10);
size_t n_eigenvectors_fine_used = ptree.get("geneo.n_eigenvectors_fine_used",(size_t)5);
size_t n_eigenvectors_coarse_computed = ptree.get("geneo.n_eigenvectors_coarse_computed",(size_t)10);
size_t n_eigenvectors_coarse_used = ptree.get("geneo.n_eigenvectors_coarse_used",(size_t)5);
size_t mlgeneo_verbose = ptree.get("geneo.verbose",(size_t)0);
double abs_zero_ker = ptree.get("geneo.abs_zero_ker",(double)1e-10);
double regularization_ker = ptree.get("geneo.regularization_ker",(double)0.0);
double eigenvalue_fine_threshold = ptree.get("geneo.eigenvalue_fine_threshold",(double)-1.0);
double eigenvalue_coarse_threshold = ptree.get("geneo.eigenvalue_coarse_threshold",(double)-1.0);
double arpack_tolerance = ptree.get("geneo.arpack_tolerance",(double)0.0);
bool merge_disconnected = ptree.get("geneo.merge_disconnected",(bool)false);
std::string solvertype = ptree.get("solver.type","direct");
// set up ISTL solver
int itmax = 100000;
using Operator = Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV>;
Operator op(istlA);
// Dune::SeqILU<ISTLM,ISTLV,ISTLV> prec(istlA,0,1.0,false);
Dune::Richardson<ISTLV,ISTLV> prec(1.0);
Dune::CGSolver<ISTLV> solver(op,prec,1e-8,itmax,2);
//Dune::RestartedGMResSolver<ISTLV> solver(op,prec,1e-8,50,itmax,2);
Dune::InverseOperatorResult stat;
Vector v(vgfs);
ISTLV& istlv = Dune::PDELab::Backend::native(v);
v = 1.0;
set_constrained_dofs(cc,0.0,v); // dirichlet dofs are zero, others are 1, so we can use it as a mask!
solver.apply(istlv,istlr,stat);
istldisplacement -= istlv;
if (solvertype=="geneo")
{
timer.reset();
MultiLevelGenEO<ISTLM,ISTLV> precx(Dune::MPIHelper::getLocalCommunicator(),
pistlA,
pistldirichlet_mask,
pvolume_snippet_matrices,
pskeleton_snippet_matrices,
pdddm,
subdomains,
pum,
fineGEVPrhs,
coarseGEVPrhs,
coarse_overlap,
coarseeigensolver,
arpack_tolerance,
n_eigenvectors_fine_computed,n_eigenvectors_fine_used,
n_eigenvectors_coarse_computed,n_eigenvectors_coarse_used,
eigenvalue_fine_threshold,
eigenvalue_coarse_threshold,
abs_zero_ker,
regularization_ker,
merge_disconnected,
cycle,
mlgeneo_verbose);
auto setup_time = timer.elapsed();
std::cout << "SETUP TIME=" << setup_time << std::endl;
// set up ISTL solver
size_t itmax = ptree.get("solver.itmax",(size_t)100);
double red = ptree.get("solver.reduction",(double)1e-6);
using Operator = Dune::MatrixAdapter<ISTLM,ISTLV,ISTLV>;
Operator op(istlA);
//Dune::SeqILU<ISTLM,ISTLV,ISTLV> prec(istlA,0,1.0,false);
Dune::Richardson<ISTLV,ISTLV> prec(1.0);
Dune::CGSolver<ISTLV> solver(op,precx,red,itmax,2);
//Dune::RestartedGMResSolver<ISTLV> solver(op,precx,red,50,itmax,2);
Dune::InverseOperatorResult stat;
Vector v(vgfs);
ISTLV& istlv = Dune::PDELab::Backend::native(v);
v = 1.0;
set_constrained_dofs(cc,0.0,v); // dirichlet dofs are zero, others are 1, so we can use it as a mask!
solver.apply(istlv,istlr,stat);
istldisplacement -= istlv;
}
//==================
// End Geneo stuff
// The old solver
//==================
// The old solver
// slp.apply();
if (solvertype=="direct")
{
using LinearSolver = Dune::PDELab::ISTLBackend_SEQ_UMFPack;
using StationaryLinearProblemSolver = Dune::PDELab::StationaryLinearProblemSolver<GridOperator, LinearSolver, Vector>;
LinearSolver ls(0);
StationaryLinearProblemSolver slp(go, ls, displacement, 1e-12);
slp.apply();
}
// Visualization
Dune::SubsamplingVTKWriter vtkwriter(gv, Dune::RefinementIntervals(2));
//Dune::VTKWriter vtkwriter(gv);
// Add the displacement field visualization
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, vgfs, displacement , Dune::PDELab::vtk::DefaultFunctionNameGenerator("Displacement"));
......@@ -223,6 +295,12 @@ int main(int argc, char** argv)
VonMisesStressGridFunction<Vector, 2> stress(displacement, material);
Dune::PDELab::interpolate(stress, p1gfs, stress_container);
Dune::PDELab::addSolutionToVTKWriter(vtkwriter, p1gfs, stress_container, Dune::PDELab::vtk::DefaultFunctionNameGenerator("von-Mises stress"));
// vtkwriter.addCellData(pdd->partition,"partition 0");
// vtkwriter.addCellData(pdd->vizpartition,"vizpartition 0");
// std::vector<unsigned> k(gv.size(0));
// for (size_t i=0; i<gv.size(0); ++i)
// k[i] = pdd->overlappingsubdomains[i].size();
// vtkwriter.addCellData(k,"k_0");
// Write visualization to file
vtkwriter.write("output", Dune::VTK::appendedraw);
......
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