Commit ff171897 authored by m.tavelli's avatar m.tavelli

Add the first version of FOCCZ4, just container, do nothing up to now

parent 89c7a876
This diff is collapsed.
This diff is collapsed.
#ifndef __FOCCZ4Solver_ADERDG_CLASS_HEADER__
#define __FOCCZ4Solver_ADERDG_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 "AbstractFOCCZ4Solver_ADERDG.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
namespace FOCCZ4{
class FOCCZ4Solver_ADERDG;
}
class FOCCZ4::FOCCZ4Solver_ADERDG : public FOCCZ4::AbstractFOCCZ4Solver_ADERDG {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
FOCCZ4Solver_ADERDG(
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
/**
* Compute the flux tensor.
*
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
*
* @param[inout] F flux at that point;
* range[outer->inner]: [0,nDim-1]x[0,nVar-1],
* already allocated.
*/
void flux(const double* const Q,double** const F) final override;
/* 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 gradQ,double* const BgradQ) final override;
void static fusedSource(const double* const restrict Q, const double* const restrict gradQ, double* const restrict S);
bool isPhysicallyAdmissible(
const double* const solution,
const double* const observablesMin,const double* const observablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
const double t) const override;
/* pointSource() function not included, as requested in the specification file */
/* multiplyMaterialParameterMatrix() not included, as requested in the specification file */
};
#endif // __FOCCZ4Solver_ADERDG_CLASS_HEADER__
\ No newline at end of file
This diff is collapsed.
#ifndef __FOCCZ4Solver_FV_CLASS_HEADER__
#define __FOCCZ4Solver_FV_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 "AbstractFOCCZ4Solver_FV.h"
#include "exahype/parser/ParserView.h"
/**
* We use Peano's logging
*/
#include "tarch/logging/Log.h"
namespace FOCCZ4{
class FOCCZ4Solver_FV;
}
class FOCCZ4::FOCCZ4Solver_FV : public FOCCZ4::AbstractFOCCZ4Solver_FV {
private:
/**
* Log device
*/
static tarch::logging::Log _log;
public:
FOCCZ4Solver_FV(
const double maximumMeshSize,
const exahype::solvers::Solver::TimeStepping timeStepping
);
/**
* 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 adjustSolution(const double* const x,const double t,const double dt, double* const Q) 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) 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 boundary; 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 vector of state variables (plus parameters) from inside of the domain
* range: [0,nVar+nPar-1], already allocated.
* @param[inout] QOut vector of state variables (plus parameters) from outside of the domain.
* range: [0,nVar+nPar-1], already allocated.
*
* @note The argument QOut is initially set as QIn mirrored at the boundary face. This
* makes it easier to implement certain boundary conditions where a velocity, e.g.,
* needs to change its sign compared to the inside state.
*/
void boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int direction,const double* const QIn,double* const QOut) override;
/**
* Compute the flux tensor.
*
* @param[in] Q vector of state variables (plus parameters);
* range: [0,nVar+nPar-1], already allocated.
*
* @param[inout] F flux at that point;
* range[outer->inner]: [0,nDim-1]x[0,nVar-1],
* already allocated.
*/
void flux(const double* const Q,double** const F) override;
/* 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] 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.
*
* @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 gradQ,double* const BgradQ) override;
/* pointSource() function not included, as requested in the specification file */
};
#endif // __FOCCZ4Solver_FV_CLASS_HEADER__
\ No newline at end of file
! DIM Initial Data
#ifndef HEADER_INITALDATA
#define HEADER_INITALDATA
!#include "MainVariables.f90"
!#include "expintegrator_type.f90"
RECURSIVE SUBROUTINE InitParameters(STRLEN,PARSETUP)
USE MainVariables, ONLY: nVar , nDim, EQN, ICType
IMPLICIT NONE
integer :: STRLEN
character(len=STRLEN) :: PARSETUP
ICType=trim(parsetup)
print *, "****************************************************************"
print *, 'Chosen setup: ',ICType
print *, "****************************************************************"
!stop
select case(ICType)
case('NLOPRUPTURE')
end select
END SUBROUTINE InitParameters
RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
USE, INTRINSIC :: ISO_C_BINDING
USE MainVariables, ONLY : nVar, nDim, EQN, ICType
IMPLICIT NONE
! Argument list
REAL, INTENT(IN) :: xGP(3), tGp !
REAL, INTENT(OUT) :: Q(nVar) !
REAL :: up(nVar)
select case(ICType)
case('NLOPRUPTURE')
CASE DEFAULT
PRINT *, 'NOT IMPLEMENTED'
STOP
END SELECT
CALL PDEPrim2Cons(Q,up)
END SUBROUTINE InitialData
RECURSIVE SUBROUTINE PDElimitervalue(limiter_value,xx,numberOfObservables, observablesMin, observablesMax)
USE, INTRINSIC :: ISO_C_BINDING
USE MainVariables, ONLY : nVar, nDim
IMPLICIT NONE
! Argument list
REAL, INTENT(IN) :: xx(3) !
INTEGER, INTENT(IN) :: numberOfObservables
INTEGER, INTENT(OUT) :: limiter_value !
REAL, INTENT(IN) :: observablesMin(numberOfObservables), observablesMax(numberOfObservables)
real :: rr
limiter_value=0
IF((observablesMax(1)<0.999 .and. observablesMin(1)>0.001) .or. observablesMax(1)>1.001) THEN ! was -0.001 .or. levelmin<0.001
limiter_value = 1
ENDIF
IF(observablesMax(2)>1.e-5 .and. observablesMax(1)>1.e-4) THEN
limiter_value = 1
ENDIF
IF(observablesMax(3)>0.5 ) THEN
limiter_value = 1
ENDIF
!limiter_value = 1
END SUBROUTINE PDElimitervalue
#endif
#ifndef __INITIAL_DATA_ADAPTER_CPP_FORTRAN_MHD__
#define __INITIAL_DATA_ADAPTER_CPP_FORTRAN_MHD__
extern "C" {
// FORTRAN functions called by C
//void initialdata_(const double* x, const double* const t, double* Q, int* md, double * cms, const int* const order);
void initialdata_(const double* x, const double* const t, double* Q);
void initparameters_(const int* length, const char* parsetup);
// only initialdata_ is used, no more.
void pdelimitervalue_(int* limiter_value, const double* xx,const int* const numberOfObservables, const double* const observablesMin,const double* const observablesMax);
}/* extern "C" */
#endif /* __INITIAL_DATA_ADAPTER_CPP_FORTRAN_MHD__ */
#ifndef HEADER_MAINVARIABLES
#define HEADER_MAINVARIABLES
! DIM Parameters
MODULE MainVariables
IMPLICIT NONE
PUBLIC
! This typesDef.f90 is a relict still from the old Fortran interface in Exahype.
! However, the following two variables are needed in the Fortran code. They
! should provided by the glue code generator in an appropriate way.
! If you modify the SRHD.exahype, please do the following mapping by hand:
!
! solver ADER-DG SRHDSolver / unknowns -> goes to -> nVar
! computational-domain / dimension -> goes to -> nDim
!
! Here, we obtain DIMENSIONS from the C++ preprocessor
#if defined(Dim3)
INTEGER, PARAMETER :: nDim = 3 ! The number of space dimensions
#elif defined(Dim2)
INTEGER, PARAMETER :: nDim = 2 ! The number of space dimensions
#endif
INTEGER, PARAMETER :: nAux = 0
INTEGER, PARAMETER :: nVar = 96 ! The number of variables of the PDE system
INTEGER, PARAMETER :: nLin = 0
INTEGER, PARAMETER :: nParam=0
INTEGER, PARAMETER :: d=3
!CHARACTER(LEN=20), PARAMETER :: ICType='NLOPRUPTURE'
CHARACTER(LEN=20) :: ICType
TYPE tEquations
REAL(8) :: gamma, Pi, c0, g = 9.81, friction = 1.0
REAL(8) :: CCZ4k1, CCZ4k2, CCZ4k3, CCZ4eta, CCZ4itau, CCZ4f, CCZ4g, CCZ4xi, CCZ4e, CCZ4c, CCZ4mu, CCZ4ds, CCZ4sk, CCZ4bs
REAL(8) :: CCZ4GLMc0 = 0.5, CCZ4GLMc = 0.75, CCZ4GLMd = 0.75, CCZ4GLMepsD = 1e-2, CCZ4GLMepsA = 1e-2, CCZ4GLMepsP = 1e-2, cs, alpha, beta, lambda, cv, rho0, p0, tau1, tau2, mu, kappa ,tau
INTEGER :: CCZ4LapseType, EinsteinAutoAux = 0, ReferenceDepth = 1.0
REAL(8) :: DivCleaning_a = 1.0
END TYPE tEquations
TYPE(tEquations) :: EQN
! 3-point Gaussian quadrature
REAL, PARAMETER :: sGP3(3) = (/ 0.5-sqrt(15.)/10.,0.5,0.5+sqrt(15.)/10. /)
REAL, PARAMETER :: wGP3(3) = (/ 5./18., 8./18., 5./18. /)
INTEGER, PARAMETER :: nGP3 = 3
END MODULE MainVariables
#endif
This diff is collapsed.
#ifndef __EXAHYPE_USER_PDE__
#define __EXAHYPE_USER_PDE__
// Fortran functions:
extern "C" {
void minimumtreedepth_(int* depth);
void hastoadjustsolution_(double* t, bool* refine);
void adjustedsolutionvalues_(const double* const x,const double* w,const double* t,const double* dt,double* Q);
void pdeflux_(double* Fx, double* Fy, double* Fz, const double* const Q);
void pdesource_(double* S, const double* const Q);
void pdencp_(double* BgradQ, const double* const Q, const double* const gradQ);
void pdefusedsrcncp_(double* S, const double* const Q, const double* const gradQ)
void pdeeigenvalues_(double* lambda, const double* const Q, double* nv);
void pdevarname_(char* MyNameOUT, int* ind);
void registerinitialdata_(const char* const id_name, int* id_name_len);
//void getnumericalsolution_(double* V,double* Q);
//void getexactsolution_(double* V,double* pos,double* timeStamp);
//void inittecplot_(const int* N_in,const int* M_in);
void pdeauxvar_(double* aux, const double* const Q,double* x, const double* const time);
void pdecritialstress_(double* CS, const double* const Q);
//void hllemfluxfv_(double* FL, double* FR, const double* const QL, const double* const QR, const double* const QavL, const double* const QavR, const int* normalNonZeroIndex);
void hllemfluxfv_(double* FL, double* FR, const double* const QL, const double* const QR, const int* normalNonZeroIndex);
void hllemriemannsolver_(const int* basisSize, const int* normalNonZeroIndex, double* FL, double* FR, const double* const QL, const double* const QR, const double* const QavL, const double* const QavR);
}/* extern "C" */
#endif /* __EXAHYPE_USER_PDE__ */
{
"project_name": "FOCCZ4",
"compiler_flags": "-DCCZ4EINSTEIN -DGLMROT",
"paths": {
"peano_kernel_path": "./Peano",
"exahype_path": "./ExaHyPE",
"output_directory": "./ApplicationExamples/FOCCZ4/FOCCZ4",
"log_file": "whatever.log"
},
"architecture": "skx",
"computational_domain": {
"dimension": 2,
"end_time": 3.0,
"offset": [
-20.0,
-20.0
],
"width": [
40.0,
40.0
]
},
"shared_memory": {
"cores": 10,
"properties_file": "sharedmemory.properties",
"autotuning_strategy": "dummy",
"background_job_consumers": 9
},
"distributed_memory": {
"timeout": 6000,
"load_balancing_type": "static",
"buffer_size": 6400,
"load_balancing_strategy": "hotspot",
"node_pool_strategy": "fair",
"ranks_per_node": 10
},
"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_update_as_background_thread": true,
"spawn_amr_background_threads": true,
"disable_vertex_exchange_in_time_steps": true,
"time_step_batch_factor": 0.0,
"disable_metadata_exchange_in_batched_time_steps": false,
"double_compression": 0.0,
"spawn_double_compression_as_background_thread": true
},
"solvers": [
{
"type": "Limiting-ADER-DG",
"name": "FOCCZ4Solver",
"order": 3,
"maximum_mesh_size": 1.0,
"maximum_mesh_depth": 2,
"time_stepping": "global",
"aderdg_kernel": {
"language": "C",
"nonlinear": true,
"terms": [
"flux",
"ncp",
"source"
],
"space_time_predictor": {},
"optimised_terms": [
"fusedsource"
],
"optimised_kernel_debugging": [],
"implementation": "optimised"
},
"point_sources": 0,
"limiter": {
"dmp_observables": 3,
"dmp_relaxation_parameter": 1e+3,
"dmp_difference_scaling": 1e+4,
"patch_size": "max",
"implementation": "generic"
},
"fv_kernel": {
"language": "C",
"terms": [
"flux",
"ncp",
"source"
],
"scheme": "musclhancock",
"slope_limiter" : "minmod",
"implementation": "generic"
},
"variables": [
{
"name": "MYQ",
"multiplicity": 96
}
],
"parameters": {
"reference": "NLOPRUPTUREPW"
},
"plotters": [
]
}
]
}
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