10.12., 9:00 - 11:00: Due to updates GitLab may be unavailable for some minutes between 09:00 and 11:00.

Commit d1c6f3e6 authored by Jean-Matthieu Gallard's avatar Jean-Matthieu Gallard

Merge branch 'master' into jm/eigen

parents 92f88d59 3dea4387
......@@ -147,4 +147,4 @@
"architecture": "noarch",
"compiler_flags": "",
"linker_flags": ""
}
\ No newline at end of file
}
......@@ -2,7 +2,7 @@
#include "Scenarios/Scenario.h"
#ifdef OPT_KERNELS
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/Kernels.h"
#include "kernels/NavierStokes_NavierStokesSolver_ADERDG/aderdg/Kernels.h"
#endif
NavierStokes::PDE::PDE() :
......
......@@ -64,4 +64,4 @@ NavierStokes::ScenarioConfig NavierStokes::parseConfig(
auto ns = PDE(referenceViscosity, *scenario, useGravity, useBackgroundState);
return {ns, scenarioName, std::move(scenario), amrSettings};
}
\ No newline at end of file
}
// This file was generated by the ExaHyPE toolkit.
// It will NOT be regenerated or overwritten.
// Please adapt it to your own needs.
//
// ========================
// www.exahype.eu
// ========================
#include "PoroelasticSolver.h"
#include "PoroelasticSolver_Variables.h"
#include "kernels/KernelUtils.h"
#include "peano/utils/Loop.h"
tarch::logging::Log Poroelastic::PoroelasticSolver::_log( "Poroelastic::PoroelasticSolver" );
void Poroelastic::PoroelasticSolver::init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// @todo Please implement/augment if required
initPointSourceLocations(cmdlineargs,constants);
}
void Poroelastic::PoroelasticSolver::adjustPointSolution(const double* const x,const double t,const double dt,double* const Q) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
if (tarch::la::equals(t,0.0)) {
Q[0] = 0.0;
Q[1] = 0.0;
Q[2] = 0.0;
Q[3] = 0.0;
Q[4] = 0.0;
Q[5] = 0.0;
Q[6] = 0.0;
Q[7] = 0.0;
Q[8] = 0.0;
Q[9] = 0.0;
Q[10] = 0.0;
Q[11] = 0.0;
Q[12] = 0.0;
Q[13] = 0.0;
Q[14] = 0.0;
Q[15] = 0.0;
Q[16] = 0.0;
Q[17] = 0.0;
Q[18] = 0.0;
Q[19] = 0.0;
Q[20] = 0.0;
Q[21] = 0.0;
Q[22] = 0.0;
}
}
void Poroelastic::PoroelasticSolver::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
stateOut[0] = 0.0;
stateOut[1] = 0.0;
stateOut[2] = 0.0;
stateOut[3] = 0.0;
stateOut[4] = 0.0;
stateOut[5] = 0.0;
stateOut[6] = 0.0;
stateOut[7] = 0.0;
stateOut[8] = 0.0;
stateOut[9] = 0.0;
stateOut[10] = 0.0;
stateOut[11] = 0.0;
stateOut[12] = 0.0;
fluxOut[0] = 0.0;
fluxOut[1] = 0.0;
fluxOut[2] = 0.0;
fluxOut[3] = 0.0;
fluxOut[4] = 0.0;
fluxOut[5] = 0.0;
fluxOut[6] = 0.0;
fluxOut[7] = 0.0;
fluxOut[8] = 0.0;
fluxOut[9] = 0.0;
fluxOut[10] = 0.0;
fluxOut[11] = 0.0;
fluxOut[12] = 0.0;
}
exahype::solvers::Solver::RefinementControl Poroelastic::PoroelasticSolver::refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& cellCentre,const tarch::la::Vector<DIMENSIONS,double>& cellSize,double t,const int level) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// Tip: See header file "peano/utils/Loop.h" for dimension-agnostic for loops.
// Example: Loop over all pointwise state variables (plus parameters)
//
// constexpr int sizeOfQ = NumberOfVariables+NumberOfParameters;
// dfor(i,Order+1) {
// const int iLinearised = dLinearised(i,Order+1);
// const double* const Q = luh + iLinearised * sizeOfQ; // pointwise state variables (plus parameters)
// // use Q[0], Q[1], ... Q[sizeOfQ-1]
// }
// @todo Please implement/augment if required
return exahype::solvers::Solver::RefinementControl::Keep;
}
//*****************************************************************************
//******************************** PDE ****************************************
// To use other PDE terms, specify them in the specification file, delete this
// file and its header and rerun the toolkit
//*****************************************************************************
void Poroelastic::PoroelasticSolver::eigenvalues(const double* const Q,const int direction,double* const lambda) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
lambda[0] = 1.0;
lambda[1] = 1.0;
lambda[2] = 1.0;
lambda[3] = 1.0;
lambda[4] = 1.0;
lambda[5] = 1.0;
lambda[6] = 1.0;
lambda[7] = 1.0;
lambda[8] = 1.0;
lambda[9] = 1.0;
lambda[10] = 1.0;
lambda[11] = 1.0;
lambda[12] = 1.0;
}
//You can either implement this method or modify fusedSource
void Poroelastic::PoroelasticSolver::algebraicSource(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, double *S) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
S[0] = 0.0;
S[1] = 0.0;
S[2] = 0.0;
S[3] = 0.0;
S[4] = 0.0;
S[5] = 0.0;
S[6] = 0.0;
S[7] = 0.0;
S[8] = 0.0;
S[9] = 0.0;
S[10] = 0.0;
S[11] = 0.0;
S[12] = 0.0;
}
void Poroelastic::PoroelasticSolver::nonConservativeProduct(const double* const Q,const double* const * const gradQ,double** const BgradQ){
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
BgradQ[0][0] = 0.0;
BgradQ[0][1] = 0.0;
BgradQ[0][2] = 0.0;
BgradQ[0][3] = 0.0;
BgradQ[0][4] = 0.0;
BgradQ[0][5] = 0.0;
BgradQ[0][6] = 0.0;
BgradQ[0][7] = 0.0;
BgradQ[0][8] = 0.0;
BgradQ[0][9] = 0.0;
BgradQ[0][10] = 0.0;
BgradQ[0][11] = 0.0;
BgradQ[0][12] = 0.0;
// @todo Please implement/augment if required
BgradQ[1][0] = 0.0;
BgradQ[1][1] = 0.0;
BgradQ[1][2] = 0.0;
BgradQ[1][3] = 0.0;
BgradQ[1][4] = 0.0;
BgradQ[1][5] = 0.0;
BgradQ[1][6] = 0.0;
BgradQ[1][7] = 0.0;
BgradQ[1][8] = 0.0;
BgradQ[1][9] = 0.0;
BgradQ[1][10] = 0.0;
BgradQ[1][11] = 0.0;
BgradQ[1][12] = 0.0;
// @todo Please implement/augment if required
BgradQ[2][0] = 0.0;
BgradQ[2][1] = 0.0;
BgradQ[2][2] = 0.0;
BgradQ[2][3] = 0.0;
BgradQ[2][4] = 0.0;
BgradQ[2][5] = 0.0;
BgradQ[2][6] = 0.0;
BgradQ[2][7] = 0.0;
BgradQ[2][8] = 0.0;
BgradQ[2][9] = 0.0;
BgradQ[2][10] = 0.0;
BgradQ[2][11] = 0.0;
BgradQ[2][12] = 0.0;
}
void Poroelastic::PoroelasticSolver::initPointSourceLocations(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants){
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
pointSourceLocation[0][0]=0.0;
pointSourceLocation[0][1]=0.0;
pointSourceLocation[0][2]=0.0;
}
void Poroelastic::PoroelasticSolver::pointSource(const double* const Q,const double* const x,const double t,const double dt, double* const forceVector, int n) {
// Tip: You find documentation for this method in header file "Poroelastic::PoroelasticSolver.h".
// Tip: See header file "Poroelastic::AbstractPoroelasticSolver.h" for toolkit generated compile-time
// constants such as Order, NumberOfVariables, and NumberOfParameters.
// @todo Please implement/augment if required
}
{
"project_name": "Poroelastic",
"paths": {
"peano_kernel_path": "./Peano",
"exahype_path" : "./ExaHyPE",
"output_directory" : "./ApplicationExamples/Poroelastic"
},
"architecture": "hsw",
"computational_domain": {
"dimension": 3,
"end_time": 1.0,
"offset": [ -20.5, 0.0 , -20.5 ],
"width" : [ 47.0, 47.0, 47.0 ]
},
"shared_memory": {
"cores": 28,
"properties_file": "sharedmemory.properties",
"autotuning_strategy": "dummy",
"background_job_consumers": 28,
"manual_pinning": false
},
"distributed_memory": {
"timeout": 6000,
"load_balancing_type": "static",
"buffer_size": 1600,
"load_balancing_strategy": "hotspot",
"node_pool_strategy": "fair",
"ranks_per_node": 1
},
"optimisation": {
"fuse_algorithmic_steps" : "all",
"fuse_algorithmic_steps_rerun_factor" : 0.99,
"fuse_algorithmic_steps_diffusion_factor" : 0.99,
"spawn_predictor_as_background_thread" : true,
"spawn_amr_background_threads" : true,
"disable_vertex_exchange_in_time_steps" : true,
"time_step_batch_factor" : 1.0,
"disable_metadata_exchange_in_batched_time_steps" : true,
"double_compression" : 0.0,
"spawn_double_compression_as_background_thread" : false
},
"solvers": [
{
"type": "ADER-DG",
"name": "PoroelasticSolver",
"order": 7,
"maximum_mesh_size": 3.0,
"maximum_mesh_depth": 0,
"time_stepping": "globalfixed",
"aderdg_kernel": {
"language": "C",
"nonlinear": false,
"terms": ["ncp", "source", "point_sources"],
"implementation": "generic",
"space_time_predictor": {"split_ck":false},
"optimised_terms": [],
"optimised_kernel_debugging": [],
"basis": "Lobatto"
},
"point_sources": 1,
"variables": [
{ "name": "sigma","multiplicity": 6 },
{ "name": "v" ,"multiplicity": 3 },
{ "name": "p" ,"multiplicity": 1 },
{ "name": "v_f" ,"multiplicity": 3 }
],
"material_parameters": [
{ "name": "bulk_solid" ,"multiplicity": 1 },
{ "name": "rho_solid","multiplicity": 1 },
{ "name": "lambda" ,"multiplicity": 1 },
{ "name": "mue","multiplicity": 1 },
{ "name": "porosity" ,"multiplicity": 1 },
{ "name": "permeability" ,"multiplicity": 1 },
{ "name": "tortuosity" ,"multiplicity": 1 },
{ "name": "bulk_fluid" ,"multiplicity": 1 },
{ "name": "rho_fluid" ,"multiplicity": 1 },
{ "name": "viscosity" ,"multiplicity": 1 }
],
"parameters": {
},
"plotters": [
{
"type": "vtu::Legendre::Vertices::Ascii",
"name": "ConservedWriter",
"time": 0.0,
"repeat": 0.1,
"output": "./conserved",
"variables": 23
}
]
}
]
}
\ No newline at end of file
#ifndef __PoroelasticSolver_CLASS_HEADER__
#define __PoroelasticSolver_CLASS_HEADER__
// This file was initially generated by the ExaHyPE toolkit.
// You can modify it in order to extend your solver with features.
// Whenever this file is present, a re-run of the ExaHyPE toolkit will
// not overwrite it. Delete it to get it regenerated.
//
// ========================
// www.exahype.eu
// ========================
#include <ostream>
#include "AbstractPoroelasticSolver.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
namespace Poroelastic{
class PoroelasticSolver;
}
class Poroelastic::PoroelasticSolver : public Poroelastic::AbstractPoroelasticSolver {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
PoroelasticSolver(
const double maximumMeshSize,
const int maximumMeshDepth,
const int haloCells,
const int haloBufferCells,
const int limiterBufferCells,
const int regularisedFineGridLevels,
const exahype::solvers::Solver::TimeStepping timeStepping,
const int DMPObservables
);
/**
* Initialise the solver.
*
* @param[in] cmdlineargs the command line arguments.
* @param[in] constants access to the constants specified for the solver.
*/
void init(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants) final override;
/**
* Adjust the conserved variables and parameters (together: Q) at a given time t at the (quadrature) point x.
*
* @note Please overwrite function adjustSolution(...) if you want to
* adjust the solution degrees of freedom in a cellwise manner.
*
* @param[in] x physical coordinate on the face.
* @param[in] t start of the time interval.
* @param[in] dt width of the time interval.
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
*/
void adjustPointSolution(const double* const x,const double t,const double dt,double* const Q) final override;
/**
* Compute the eigenvalues of the flux tensor per coordinate @p direction.
*
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
* @param[in] direction normal direction of the face / column of the flux vector (range: [0,nDims-1]).
* @param[inout] lambda eigenvalues as C array;
* range: [0,nVar-1], already allocated.
*/
void eigenvalues(const double* const Q,const int direction,double* const lambda) final override;
/**
* Impose boundary conditions at a point on a boundary face
* within the time interval [t,t+dt].
*
* @param[in] x physical coordinate on the face;
* range: [0,nDims-1].
* @param[in] t start of the time interval.
* @param[in] dt width of the time interval.
* @param[in] faceIndex indexing of the face (0 -- {x[0]=xmin}, 1 -- {x[1]=xmax}, 2 -- {x[1]=ymin}, 3 -- {x[2]=ymax}, and so on,
* where xmin,xmax,ymin,ymax are the bounds of the cell containing point x.
* @param[in] direction coordinate direction the face normal is pointing to.
*
* @param[in] QIn the conserved variables (plus parameters) at point x from inside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar+nPar-1], already allocated.
* @param[in] FIn the normal fluxes at point x from inside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar-1], already allocated.
*
* @param[inout] QOut the conserved variables at point x from outside of the domain
* and time-averaged (over [t,t+dt]);
* range: [0,nVar+nPar-1], already allocated.
* @param[inout] FOut the normal fluxes at point x from outside of the domain;
* range: [0,nVar-1], already allocated.
*/
void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double * const fluxIn,const double* const stateIn, const double* const gradStateIn,double *fluxOut,double* stateOut) final override;
/**
* Evaluate the refinement criterion within a cell.
*
* @note Instead of a variables array at a single quadrature point we give
* you all NumberOfVariables*(Order+1)^DIMENSIONS solution degrees of freedom.
*
* @note Use this function and ::adjustSolution to set initial conditions.
*
* @param[in] luh all state variables (plus parameters) of a cell;
* range[outer->inner]: [0,(p+1)^nDim-1]x[0,nVar+nPar-1],
* already allocated.
*
* @param[in] cellCentre cellCentre of the cell.
* @param[in] cellSize extent of the cell.
* @param[in] t start of the time interval.
* @param[in] dt width of the time interval.
*
* @return One of exahype::solvers::Solver::RefinementControl::{Erase,Keep,Refine}.
*/
exahype::solvers::Solver::RefinementControl refinementCriterion(const double* const luh,const tarch::la::Vector<DIMENSIONS,double>& cellCentre,const tarch::la::Vector<DIMENSIONS,double>& cellSize,double t,const int level) override;
//PDE
/* flux() function not included, as requested in the specification file */
/* viscousFlux() function not included, as requested in the specification file */
/**
* Compute the Algebraic Sourceterms.
*
* You may want to overwrite this with your PDE Source (algebraic RHS contributions).
* However, in all schemes we have so far, the source-type contributions are
* collected with non-conservative contributions into a fusedSource, see the
* fusedSource method. From the kernels given with ExaHyPE, only the fusedSource
* is called and there is a default implementation for the fusedSource calling
* again seperately the nonConservativeProduct function and the algebraicSource
* function.
*
* @param[in] x physical position
* @param[in] t physical time
* @param[in] Q vector of state variables (plus material
* parameters); range: [0,nVar+nPar-1], already allocated.
* @param[inout] S source term; range: [0,nVar-1], already allocated.
*/
void algebraicSource(const tarch::la::Vector<DIMENSIONS, double>& x, double t, const double *const Q, double *S) override;
/**
* Compute the nonconservative term $B(Q) \nabla Q$.
*
* This function shall return a vector BgradQ which holds the result
* of the full term. To do so, it gets the vector Q and the matrix
* gradQ which holds the derivative of Q in each spatial direction.
* Currently, the gradQ is a continous storage and users can use the
* kernels::icellSize2 class in order to compute the positions inside gradQ.
*
* @TODO: Check if the following is still right:
*
* !!! Warning: BgradQ is a vector of size NumberOfVariables if you
* use the ADER-DG kernels for nonlinear PDEs. If you use
* the kernels for linear PDEs, it is a tensor with dimensions
* Dim x NumberOfVariables.
*
* @param[in] Q vector of state variables (plus material
* parameters); range: [0,nVar+nPar-1], already allocated.
* @param[in] gradQ the gradients of the vector of unknowns, stored
* in a linearized array. (range: [0,dim*(nVar+nPar)-1].
* @param[inout] BgradQ the nonconservative product (extends nVar),
* already allocated.
*/
void nonConservativeProduct(const double* const Q,const double* const * const gradQ,double** const BgradQ);
/**
* Initialise the locations of point sources
* Point source n will be located at pointSourceLocation[n][:]
*
* @param[in] cmdlineargs the command line arguments.
* @param[in] constants access to the constants specified for the solver.
*/
void initPointSourceLocations(const std::vector<std::string>& cmdlineargs,const exahype::parser::ParserView& constants);
/**
* Compute the contribution of a point source.
*
* This function should return the force vector of a point source n
* located at pointSourceLocation[n][:] as defined in
* initPointSourceLocations()
*
* @param[in] Q vector of state variables (plus material
* parameters), range: [0,nVar+nPar-1], already allocated.
* @param[in] x the position associated with Q
* @param[in] t the time stamp associated with Q
* @param[in] dt the current time step size
* @param[out] forceVector the force vector of the point source (size: ?)
*
* @todo LR: Is x really needed?
* @todo LR: Is dt really needed?
* @todo LR: specify size of forceVector
*/
void pointSource(const double* const Q,const double* const x,const double t,const double dt, double* const forceVector,int n) override;
/* multiplyMaterialParameterMatrix() not included, as requested in the specification file */
};
#endif // __PoroelasticSolver_CLASS_HEADER__
......@@ -216,17 +216,6 @@ exahype::solvers::ADERDGSolver::ADERDGSolver(
_profiler->registerTag(tag);
}
#ifdef Parallel
_invalidExtrapolatedPredictor.resize(getBndFaceSize());
_invalidFluctuations.resize(getBndFluxSize());
std::fill_n(_invalidExtrapolatedPredictor.data(),_invalidExtrapolatedPredictor.size(),-1);
std::fill_n(_invalidFluctuations.data(),_invalidFluctuations.size(),-1);
_receivedExtrapolatedPredictor.resize(getBndFaceSize());
_receivedFluctuations.resize(getBndFluxSize());
_receivedUpdate.reserve(getUpdateSize());
#endif
}
int exahype::solvers::ADERDGSolver::getUnknownsPerFace() const {
......@@ -500,6 +489,19 @@ void exahype::solvers::ADERDGSolver::initSolver(
resetGlobalObservables(_nextGlobalObservables.data());
init(cmdlineargs,parserView); // call user define initalisiation
#ifdef Parallel
_invalidExtrapolatedPredictor.resize(getBndFaceSize());
_invalidFluctuations.resize(getBndFluxSize());
std::fill_n(_invalidExtrapolatedPredictor.data(),_invalidExtrapolatedPredictor.size(),-1);
std::fill_n(_invalidFluctuations.data(),_invalidFluctuations.size(),-1);
_receivedExtrapolatedPredictor.resize(getBndFaceSize());
_receivedFluctuations.resize(getBndFluxSize());
_receivedUpdate.reserve(getUpdateSize());
#endif
}
bool exahype::solvers::ADERDGSolver::isPerformingPrediction(
......@@ -2009,7 +2011,7 @@ void exahype::solvers::ADERDGSolver::mergeWithNeighbourData(
DataHeap::getInstance().receiveData( // TODO const-correct peano
const_cast<double*>(_receivedExtrapolatedPredictor.data()),dataPerFace,
fromRank, barycentre, level, peano::heap::MessageType::NeighbourCommunication);
solveRiemannProblemAtInterface(
cellDescription, face,