The TS-GitLab will have to shut down in the near future — please plan migrating your projects to GitLab.com or GitHub. Contact @yunus for more information.

Commit 733c4678 authored by Santiago Ospina's avatar Santiago Ospina
Browse files

[Flux Reconstruction] Add test to check normal jumps

parent 847f1cd5
......@@ -76,6 +76,20 @@ add_custom_target(test_mass_conservation
COMMAND ctest --output-on-failure --tests-regex ^.+mass-conservation.+$
)
dorie_add_metaini_test(
SOURCE test-conforming-flux-jumps.cc
BASENAME conforming-flux-jumps
CREATED_TARGETS conforming_flux_jumps
METAINI conforming_flux_jumps.mini.in
SCRIPT dune_execute.py
)
target_link_libraries(${conforming_flux_jumps} dorie-impl spdlog)
add_custom_target(test_conforming_flux_jumps
COMMAND ctest --output-on-failure --tests-regex ^.+conforming-flux-jumps.+$
)
# dune excludes test targets from 'make all'; undo that here where applicable
set_property(TARGET dorie PROPERTY EXCLUDE_FROM_ALL 0)
set_property(TARGET dorie-rfg PROPERTY EXCLUDE_FROM_ALL 0)
include ${CMAKE_BINARY_DIR}/doc/default_files/config.ini
_asset_path = "${CMAKE_CURRENT_LIST_DIR}"
_order = 1, 2, 3 | expand order
_dim = 2, 3 | expand dim
_entity_type = _simplex, _cube | expand grid
__name = conforming_flux_jumps_{_dim}D_{_order}{_entity_type}
({_order} == 3 and {_entity_type} == _cube) and {_dim} == 3 | exclude
_limit_conforming_flux_jumps_2D_1_cube = 1E-20
_limit_conforming_flux_jumps_2D_1_simplex = 1E-20
_limit_conforming_flux_jumps_2D_2_cube = 1E-20
_limit_conforming_flux_jumps_2D_2_simplex = 1E-20
_limit_conforming_flux_jumps_2D_3_cube = 1E-17
_limit_conforming_flux_jumps_2D_3_simplex = 1E-20
_limit_conforming_flux_jumps_3D_1_cube = 1E-15
_limit_conforming_flux_jumps_3D_1_simplex = 1E-15
_limit_conforming_flux_jumps_3D_2_cube = 1E-15
# wait... what?
_limit_conforming_flux_jumps_3D_2_simplex = 1E-4
_limit_conforming_flux_jumps_3D_3_cube = 1E-15
_limit_conforming_flux_jumps_3D_3_simplex = 1E-15
[adaptivity]
policy = none
[grid]
dimensions = {_dim}
gridType = gmsh, rectangular | expand grid
initialLevel = 0
cells = 20 20, 10 10 10 | expand dim
extensions = 1 1, 1 1 1 | expand dim
gridFile = {_asset_path}/meshes/square.msh, {_asset_path}/meshes/3dcube.msh | expand dim
[richards]
check.limitConformingFluxJumps = {_limit_{__name}}
output.fileName = {__name}
output.outputPath = {__name}
output.vertexData = true
boundary.file = "{_asset_path}/bcs/infiltration_2d.dat", "{_asset_path}/bcs/infiltration_3d.dat" | expand dim
parameters.file = "{_asset_path}/param/param.yml"
time.end = 2E2
time.maxTimestep = 2E2
time.startTimestep = 1E2
# Increase jumps in dG solution!
numerics.penaltyFactor = 1E4
numerics.FEorder = {_order}
\ No newline at end of file
// -*- tab-width: 2; indent-tabs-mode: nil -*-
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <string>
#include <iostream>
#include <cstdio>
#include <exception>
#include <unistd.h>
#include <sys/stat.h>
#include <sys/types.h>
#include <dune/common/exceptions.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/dorie/common/setup.hh>
#include <dune/dorie/common/logging.hh>
#include <dune/dorie/common/grid_creator.hh>
#include <dune/dorie/model/richards/richards.hh>
#include "test-conforming-flux-jumps.hh"
//===============================================================
// Main program with grid setup
//===============================================================
template<typename Traits>
using Sim = Dune::Dorie::TestConformingFluxJumps<Traits>;
template<int dim, int order>
using Simplex = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::UGGrid<dim>,
Dune::GeometryType::BasicType::simplex>,order>;
template<int dim, int order>
using Cube = Dune::Dorie::RichardsSimulationTraits<Dune::Dorie::BaseTraits<Dune::YaspGrid<dim>,
Dune::GeometryType::BasicType::cube>,order>;
/// Main Program Function: Read config file, build grid and call Richards Solver
int main(int argc, char** argv)
{
try{
Dune::Timer timer;
// Initialize all the stuff!
const std::string grt = "Starting DORiE";
auto [inifile, log, helper] = Dune::Dorie::Setup::init(argc, argv, grt);
// halt process for debugger
if (inifile.get<bool>("misc.debugMode")) {
Dune::Dorie::Setup::debug_hook();
}
// setup richards configuration
Dune::ParameterTree richards_config = Dune::Dorie::Setup::prep_ini_for_richards(inifile);
// Read necessary variables
const std::string gtype = inifile.get<std::string>("grid.gridType");
const int dim = inifile.get<int>("grid.dimensions");
const std::string adapt_policy_str = inifile.get<std::string>("adaptivity.policy");
const int FEorder = richards_config.get<int>("numerics.FEorder");
const std::string outputPath = richards_config.get<std::string>("output.outputPath");
Dune::Dorie::AdaptivityPolicy adapt_policy = Dune::Dorie::AdaptivityPolicy::None;
if (adapt_policy_str == "waterFlux")
adapt_policy = Dune::Dorie::AdaptivityPolicy::WaterFlux;
else if (adapt_policy_str != "none")
DUNE_THROW(Dune::NotImplemented,"The adaptivity policy " <<
adapt_policy_str << " is not implemented!");
// Attempt to create output directory
log->info("Creating output directory: {}", outputPath);
mkdir(outputPath.c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
int result = access(outputPath.c_str(), W_OK);
if (result != 0)
DUNE_THROW(Dune::IOError,"Output folder " << outputPath << " not writable");
if (dim==2)
{
if (gtype == "gmsh"){
Dune::Dorie::GridCreator<Dune::UGGrid<2>> grid_creator(inifile, helper);
switch(FEorder){
case 1:{
Sim<Simplex<2,1>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<Simplex<2,2>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{ // No flux reconstruction available
Sim<Simplex<2,3>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
if (adapt_policy == Dune::Dorie::AdaptivityPolicy::None) {
Dune::Dorie::GridCreator<Dune::YaspGrid<2>> grid_creator(inifile, helper);
switch(FEorder){
case 1:{
Sim<Cube<2,1>> sim(richards_config, grid_creator, helper);
sim.run();
break;
}
case 2:{
Sim<Cube<2,2>> sim(richards_config, grid_creator, helper);
sim.run();
break;
}
case 3:{
Sim<Cube<2,3>> sim(richards_config, grid_creator, helper);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
} else
DUNE_THROW(Dune::NotImplemented,"Cube adaptive grids not available!");
}
else
DUNE_THROW(Dune::NotImplemented,"Grid Type not supported!");
}
else if (dim==3)
{
if (gtype == "gmsh"){
Dune::Dorie::GridCreator<Dune::UGGrid<3>> grid_creator(inifile, helper);
switch(FEorder){
case 1:{
Sim<Simplex<3,1>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 2:{
Sim<Simplex<3,2>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
case 3:{
Sim<Simplex<3,3>> sim(richards_config, grid_creator, helper);
sim.set_policy(adapt_policy);
sim.run();
break;
}
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
}
else if (gtype == "rectangular"){
if(adapt_policy == Dune::Dorie::AdaptivityPolicy::None){
Dune::Dorie::GridCreator<Dune::YaspGrid<3>> grid_creator(inifile, helper);
switch(FEorder){
case 1:{
Sim<Cube<3,1>> sim(richards_config, grid_creator, helper);
sim.run();
break;
}
case 2:{
Sim<Cube<3,2>> sim(richards_config, grid_creator, helper);
sim.run();
break;
}
// case 3:{ // No flux reconstruction available
// Sim<Cube<3,3>> sim(richards_config, grid_creator, helper);
// sim.run();
// break;
// }
default:
DUNE_THROW(Dune::NotImplemented,"Finite Element Order (numerics.FEorder) not supported!");
}
} else
DUNE_THROW(Dune::NotImplemented,"Cube adaptive grids not available!");
}
else
DUNE_THROW(Dune::NotImplemented,"Grid Type not supported!");
}
// grid_dim != 2,3
else{
DUNE_THROW(Dune::NotImplemented,"Number of dimensions (grid.dimensions) not supported!");
}
log->info("DORiE finished after {:.2e}s :)",
timer.elapsed());
return 0;
}
catch (Dune::Exception &exc){
auto log = Dune::Dorie::create_base_logger();
log->critical("Aborting DORiE after exception: {}", exc.what());
return 1;
}
catch (std::exception& exc) {
auto log = Dune::Dorie::create_base_logger();
log->critical("Aborting DORiE after exception: {}", exc.what());
return 1;
}
catch (...){
auto log = Dune::Dorie::create_base_logger();
log->error("Unknown error occurred");
log->critical("Aborting DORiE after uncaught exception");
return 1;
}
}
#ifndef DUNE_DORIE_TEST_CONFIRMING_FLUX_JUMPS_HH
#define DUNE_DORIE_TEST_CONFIRMING_FLUX_JUMPS_HH
#include <dune/dorie/model/richards/richards.cc>
namespace Dune {
namespace Dorie {
template<class Traits>
class TestConformingFluxJumps : public RichardsSimulation<Traits>
{
public:
using RichardsSimulation<Traits>::RichardsSimulation;
virtual void step() override
{
double tol = this->inifile.template get<double>("check.limitConformingFluxJumps");
RichardsSimulation<Traits>::step();
if (this->flux_rec_gf) {
if (not this->flux_rec_gf->check_jumps(tol)) {
DUNE_THROW(MathError,"Fluxes are not conforming!");
}
} else {
DUNE_THROW(InvalidStateException,"Flux reconstruction pointer is in invalid state!");
}
}
};
} // namespace Dorie
} // namespace Dune
#endif // DUNE_DORIE_TEST_CONFIRMING_FLUX_JUMPS_HH
......@@ -285,10 +285,6 @@ public:
acc += deviation;
acc_square += deviation * deviation;
if (flux_rec_gf->check_jumps(1e-10))
DUNE_THROW(Dune::MathError,
"Jumps in flux reconstruction are not admissible!");
std::cout << "wc_old: " << wc_old << std::endl;
std::cout << "wc_new: " << wc_new << std::endl;
std::cout << "integrated flux: " << flux_int << std::endl;
......
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