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

Commit f9ac52b5 authored by Philipp Samfaß's avatar Philipp Samfaß

got Euler Flow + GRMHD working again (thanks to Anne)

parent c8f0e418
#include "ComputeGlobalIntegrals.h"
#include "Primitives.h"
#include "InitialData.h"
Euler::ComputeGlobalIntegrals::ComputeGlobalIntegrals(MyEulerSolver& solver) :
conserved("output/conserved-"),
primitives("output/primitive-"),
errors("output/error-"),
statistics("output/volume.asc")
{
conserved.add(0, "dens");
conserved.add(1, "sconx");
conserved.add(2, "scony");
conserved.add(3, "sconz");
conserved.add(4, "ener");
primitives.add(0, "rho"); // V[0]=Q[0]
primitives.add(1, "velx");
primitives.add(2, "vely");
primitives.add(3, "velz");
primitives.add(4, "press");
errors.add(0, "rho");
errors.add(1, "velx");
errors.add(2, "vely");
errors.add(3, "velz");
errors.add(4, "press");
// here, we plot all variables.
assert( conserved.size() <= NumberOfVariables );
assert( primitives.size() <= NumberOfVariables );
assert( errors.size() <= NumberOfVariables );
}
Euler::ComputeGlobalIntegrals::~ComputeGlobalIntegrals() {
}
void Euler::ComputeGlobalIntegrals::startPlotting(double time) {
conserved.startRow(time);
primitives.startRow(time);
errors.startRow(time);
statistics.startRow(time);
}
void Euler::ComputeGlobalIntegrals::finishPlotting() {
conserved.finishRow();
primitives.finishRow();
errors.finishRow();
statistics.finishRow();
}
void Euler::ComputeGlobalIntegrals::mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp
) {
// make sure this plotter has no output associated
assertion( outputQuantities == nullptr );
const double NumberOfLagrangePointsPerAxis = Euler::MyEulerSolver::Order + 1;
//const double NumberOfUnknownsPerGridPoint = NumberOfVariables;
// volume form for integration
double scaling = tarch::la::volume(sizeOfPatch* (1.0/NumberOfLagrangePointsPerAxis));
statistics.addValue(scaling, 1);
// NOTE: This is *not* the correct scaling for the ADERDG scheme.
// Instead, you should use the GlobalIntegralsLegendre plotter.
// reduce the conserved quantities
conserved.addValue(Q, scaling);
// reduce the primitive quantities
double V[NumberOfVariables];
cons2prim(V, Q);
primitives.addValue(V, scaling);
// now do the convergence test, as we have exact initial data
double ExactCons[NumberOfVariables];
double ExactPrim[NumberOfVariables];
const double *xpos = x.data();
idfunc(xpos, ExactCons, timeStamp);
cons2prim(ExactPrim, ExactCons);
double localError[NumberOfVariables];
for(int i=0; i<NumberOfVariables; i++) {
localError[i] = std::abs(V[i] - ExactPrim[i]);
}
errors.addValue(localError, scaling);
}
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
#ifndef __ComputeGlobalIntegrals_CLASS_HEADER__
#define __ComputeGlobalIntegrals_CLASS_HEADER__
#include "exahype/plotters/Plotter.h"
namespace Euler{
class ComputeGlobalIntegrals;
/**
* Forward declaration
*/
class MyEulerSolver;
}
#include "TimeSeriesReductions.h"
#include "MyEulerSolver.h"
#ifndef __NumberOfVariables__
static const int NumberOfVariables = Euler::MyEulerSolver::NumberOfVariables; // shortcut
#define __NumberOfVariables__
#endif
class Euler::ComputeGlobalIntegrals: public exahype::plotters::Plotter::UserOnTheFlyPostProcessing{
private:
reductions::MultipleReductionsWriter conserved;
reductions::MultipleReductionsWriter primitives;
reductions::MultipleReductionsWriter errors;
reductions::ReductionsWriter statistics;
public:
ComputeGlobalIntegrals(MyEulerSolver& solver);
virtual ~ComputeGlobalIntegrals();
void startPlotting(double time) override;
void finishPlotting() override;
void mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp) override;
};
#endif
#include "ComputeGlobalIntegralsLegendre.h"
#include "Primitives.h"
#include "InitialData.h"
#include "kernels/GaussLegendreBasis.h"
Euler::ComputeGlobalIntegralsLegendre::ComputeGlobalIntegralsLegendre(MyEulerSolver& solver) :
conserved("output/conserved-"),
primitives("output/primitive-"),
errors("output/error-"),
statistics("output/volume.asc")
{
conserved.add(0, "dens");
conserved.add(1, "sconx");
conserved.add(2, "scony");
conserved.add(3, "sconz");
conserved.add(4, "ener");
primitives.add(0, "rho"); // V[0]=Q[0]
primitives.add(1, "velx");
primitives.add(2, "vely");
primitives.add(3, "velz");
primitives.add(4, "presss");
errors.add(0, "rho");
errors.add(1, "velx");
errors.add(2, "vely");
errors.add(3, "velz");
errors.add(4, "presss");
// here, we plot all variables.
assert( conserved.size() <= NumberOfVariables );
assert( primitives.size() <= NumberOfVariables );
assert( errors.size() <= NumberOfVariables );
}
Euler::ComputeGlobalIntegralsLegendre::~ComputeGlobalIntegralsLegendre() {
}
void Euler::ComputeGlobalIntegralsLegendre::startPlotting(double time) {
conserved.startRow(time);
primitives.startRow(time);
errors.startRow(time);
statistics.startRow(time);
}
void Euler::ComputeGlobalIntegralsLegendre::finishPlotting() {
conserved.finishRow();
primitives.finishRow();
errors.finishRow();
statistics.finishRow();
}
void Euler::ComputeGlobalIntegralsLegendre::mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp
) {
// make sure this plotter has no output associated
assertion( outputQuantities == nullptr );
// volume form for integration
double scaling = tarch::la::volume(sizeOfPatch);
statistics.addValue(scaling, 1);
// Gauss-Legendre weights from pos argument
double wx = kernels::legendre::weights[Euler::MyEulerSolver::Order][pos[0]];
double wy = kernels::legendre::weights[Euler::MyEulerSolver::Order][pos[1]];
double wz = 1;
#ifdef Dim3
wz = kernels::legendre::weights[Euler::MyEulerSolver::Order][pos[2]];
#endif
scaling *= wx*wy*wz;
// reduce the conserved quantities
conserved.addValue(Q, scaling);
// reduce the primitive quantities
double V[NumberOfVariables];
cons2prim(V, Q);
primitives.addValue(V, scaling);
// now do the convergence test, as we have exact initial data
double Exact[NumberOfVariables];
const double *xpos = x.data();
idfunc(xpos, Exact, timeStamp);
double ExactPrim[NumberOfVariables];
cons2prim(ExactPrim, Exact);
// Uncomment for debugging reasons
// std::cout << "x="<<x.toString()<<"J="<<scaling << std::endl;
// std::cout << "wx="<<wx<<",wy="<<wy<<",wz="<<wz << std::endl;
double localError[NumberOfVariables];
for(int i=0; i<NumberOfVariables; i++) {
localError[i] = std::abs(V[i] - ExactPrim[i]);
}
errors.addValue(localError, scaling);
// Uncomment for debugging reasons
// std::cout << "V_ana["<<i<<"]="<<ExactPrim[i]<<",V_h["<<i<<"]="<<V[i]<< std::endl;
// std::cout << "Q_ana["<<i<<"]="<<Exact[i]<<",Q_h["<<i<<"]="<<Q[i] <<std::endl;
}
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
#ifndef __ComputeGlobalIntegralsLegendre_CLASS_HEADER__
#define __ComputeGlobalIntegralsLegendre_CLASS_HEADER__
#include "exahype/plotters/Plotter.h"
namespace Euler{
class ComputeGlobalIntegralsLegendre;
/**
* Forward declaration
*/
class MyEulerSolver;
}
/* I hope these modifications are not overwritten... */
#include "TimeSeriesReductions.h"
#include "MyEulerSolver.h"
#ifndef __NumberOfVariables__
static const int NumberOfVariables = Euler::MyEulerSolver::NumberOfVariables; // shortcut
#define __NumberOfVariables__
#endif
class Euler::ComputeGlobalIntegralsLegendre: public exahype::plotters::Plotter::UserOnTheFlyPostProcessing{
private:
reductions::MultipleReductionsWriter conserved;
reductions::MultipleReductionsWriter primitives;
reductions::MultipleReductionsWriter errors;
reductions::ReductionsWriter statistics;
public:
ComputeGlobalIntegralsLegendre(MyEulerSolver& solver);
virtual ~ComputeGlobalIntegralsLegendre();
void startPlotting(double time) override;
void finishPlotting() override;
void mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double timeStamp) override;
};
#endif
......@@ -13,12 +13,14 @@ void InitialData(const double* const x, double* Q, double t);
void Polynomial(const double* const x, double* Q, double t);
//ok
void ShuVortex2D(const double* const x, double* Q, double t);
//ok
void MovingGauss(const double* const x, double* V, double t);
void DiffusingGauss(const double* const x, double* Q, double /*t*/);
//Ok ->am spannendsten, Gitter muss fein genug sein
void SmoothedSodShockTube(const double* const x, double* Q, double /*t*/);
#endif /* __INITIAL_DATA_EULERFLOW__ */
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "PrimitivesWriter.h"
#include "Primitives.h"
Euler::PrimitivesWriter::PrimitivesWriter(MyEulerSolver& solver) {
// @todo Please insert your code here
Euler::PrimitivesWriter::PrimitivesWriter(Euler::MyEulerSolver& solver) {
// @TODO Please insert your code here.
}
Euler::PrimitivesWriter::~PrimitivesWriter() {
// @todo Please insert your code here
}
void Euler::PrimitivesWriter::startPlotting(double time) {
// @todo Please insert your code here
void Euler::PrimitivesWriter::startPlotting( double time) {
// @TODO Please insert your code here.
}
void Euler::PrimitivesWriter::finishPlotting() {
// @todo Please insert your code here
// @TODO Please insert your code here.
}
void Euler::PrimitivesWriter::mapQuantities(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& x,
const tarch::la::Vector<DIMENSIONS, int>& pos,
double* Q,
double* outputQuantities,
double* const Q,
double* const outputQuantities,
double timeStamp
) {
/**
* This is the Primitive Variable plotter
**/
//double V[Euler::MyEulerSolver::NumberOfVariables];
cons2prim(outputQuantities, Q);
//for (int i=0; i<0; i++){
// outputQuantities[i] = V[i];
//}
}
const int writtenUnknowns = 5;
for (int i=0; i<writtenUnknowns; i++){
outputQuantities[i] = Q[i];
}
}
\ No newline at end of file
......@@ -73,11 +73,11 @@
REAL, PARAMETER :: ExcisionRadius=1.5
!
REAL, PARAMETER :: aom = 0.0, Mbh = 1.0
!CHARACTER(LEN=200), PARAMETER :: ICType = TRIM('GRMHD-SphericalAccretion')
!CHARACTER(LEN=200), PARAMETER :: ICType2 = TRIM('GRMHD-SphericalAccretion')
! CHARACTER(LEN=200), PARAMETER :: ICType2 = TRIM('GRMHD-BondiAccretion')
CHARACTER(LEN=200), PARAMETER :: ICType = TRIM('GRMHDTOV')
CHARACTER(LEN=200), PARAMETER :: ICType = TRIM('GRMHD-SphericalAccretion')
CHARACTER(LEN=200), PARAMETER :: ICType2 = TRIM('GRMHD-SphericalAccretion')
! CHARACTER(LEN=200), PARAMETER :: ICType2 = TRIM('GRMHD-BondiAccretion')
! CHARACTER(LEN=200), PARAMETER :: ICType = TRIM('GRMHDTOV')
! CHARACTER(LEN=200), PARAMETER :: ICType2 = TRIM('GRMHD-SphericalAccretion')
!ICType = 'GRMHDTOV'
!
END MODULE Parameters
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
// ==============================================
// Please do not change the implementations below
// =============================---==============
#include "GRMHDSolver.h"
#include "kernels/GRMHD_GRMHDSolver_ADERDG/Kernels.h"
// Just call parent constructor
GRMHD::GRMHDSolver::GRMHDSolver(
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,
const double DMPRelaxationParameter,
const double DMPDifferenceScaling
) :
exahype::solvers::LimitingADERDGSolver::LimitingADERDGSolver(
"GRMHDSolver",
new GRMHD::GRMHDSolver_ADERDG(
maximumMeshSize,maximumMeshDepth,haloCells,haloBufferCells,limiterBufferCells,regularisedFineGridLevels,timeStepping,DMPObservables),
new GRMHD::GRMHDSolver_FV(
maximumMeshSize, timeStepping),
DMPRelaxationParameter,
DMPDifferenceScaling) {}
void GRMHD::GRMHDSolver::projectOnFVLimiterSpace(const double* const luh, double* const lim) const {
GRMHD::GRMHDSolver_ADERDG_kernels::aderdg::projectOnFVLimiterSpace(luh, lim);
}
void GRMHD::GRMHDSolver::projectOnDGSpace(const double* const lim, double* const luh) const {
GRMHD::GRMHDSolver_ADERDG_kernels::aderdg::projectOnDGSpace(lim, luh);
}
bool GRMHD::GRMHDSolver::discreteMaximumPrincipleAndMinAndMaxSearch(const double* const luh, double* const boundaryMinPerVariables, double* const boundaryMaxPerVariables) {
const bool r = GRMHD::GRMHDSolver_ADERDG_kernels::aderdg::discreteMaximumPrincipleAndMinAndMaxSearch(luh, _solver.get(), _DMPMaximumRelaxationParameter, _DMPDifferenceScaling, boundaryMinPerVariables, boundaryMaxPerVariables);
return r;
}
void GRMHD::GRMHDSolver::findCellLocalMinAndMax(const double* const luh, double* const localMinPerVariables, double* const localMaxPerVariable) {
GRMHD::GRMHDSolver_ADERDG_kernels::aderdg::findCellLocalMinAndMax(luh, _solver.get(), localMinPerVariables, localMaxPerVariable);
}
void GRMHD::GRMHDSolver::findCellLocalLimiterMinAndMax(const double* const lim, double* const localMinPerObservable, double* const localMaxPerObservable) {
GRMHD::GRMHDSolver_ADERDG_kernels::aderdg::findCellLocalLimiterMinAndMax(lim, _solver.get(), localMinPerObservable,localMaxPerObservable);
}
\ No newline at end of file
#ifndef __GRMHDSolver_CLASS_HEADER__
#define __GRMHDSolver_CLASS_HEADER__
// This file is generated by the ExaHyPE toolkit.
// Please do not modify - it will be overwritten by the next
// ExaHyPE toolkit call.
//
// ========================
// www.exahype.eu
// ========================
#include <string>
#include "exahype/solvers/LimitingADERDGSolver.h"
#include "GRMHDSolver_ADERDG.h"
#include "GRMHDSolver_FV.h"
namespace GRMHD{
class GRMHDSolver;
}
class GRMHD::GRMHDSolver: public exahype::solvers::LimitingADERDGSolver {
public:
static constexpr int NumberOfVariables = GRMHD::AbstractGRMHDSolver_ADERDG::NumberOfVariables;
static constexpr int NumberOfParameters = GRMHD::AbstractGRMHDSolver_ADERDG::NumberOfParameters;
static constexpr int Order = GRMHD::AbstractGRMHDSolver_ADERDG::Order;
static constexpr int NumberOfDMPObservables = GRMHD::AbstractGRMHDSolver_ADERDG::NumberOfDMPObservables;
static constexpr int GhostLayerWidth = GRMHD::AbstractGRMHDSolver_FV::GhostLayerWidth;
GRMHDSolver(
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,
const double DMPRelaxationParameter,
const double DMPDifferenceScaling
);
void projectOnFVLimiterSpace(const double* const luh, double* const lim) const override;
void projectOnDGSpace(const double* const lim, double* const luh) const override;
bool discreteMaximumPrincipleAndMinAndMaxSearch(const double* const luh, double* const boundaryMinPerVariables, double* const boundaryMaxPerVariables) override;
void findCellLocalMinAndMax(const double* const luh, double* const localMinPerVariables, double* const localMaxPerVariable) override;
void findCellLocalLimiterMinAndMax(const double* const lim, double* const localMinPerObservable, double* const localMaxPerObservable) override;
};
#endif // __GRMHDSolver_CLASS_HEADER__
\ No newline at end of file
......@@ -33,10 +33,10 @@ void GRMHD::GRMHDSolver_ADERDG::init(const std::vector<std::string>& cmdlineargs
// feenableexcept(FE_INVALID | FE_OVERFLOW); // Enable all floating point exceptions but FE_INEXACT
// Todo: Move this to specfile once we have working constants.
// std::string id_default = "Fortran";
// std::string bc_default = "left:exact,right:exact,top:exact,bottom:exact,front:exact,back:exact";
std::string id_default = "TOVSolver";
std::string bc_default = "left:exact,right:exact,top:exact,bottom:exact,front:exact,back:exact";
std::string id_default = "Fortran";
std::string bc_default = "left:exact,right:exact,top:exact,bottom:exact,front:exact,back:exact";
//std::string id_default = "TOVSolver";
//std::string bc_default = "left:exact,right:exact,top:exact,bottom:exact,front:exact,back:exact";
// alternatives:
//std::string id_default = "RNSID";
......
......@@ -13,7 +13,7 @@ bool prepare_id(std::string idname) {
if(idname == "Fortran") {
id = new fortranid();
return true;
} else if(idname == "PizzaTOV") {
}/* else if(idname == "PizzaTOV") {
id = new pizzatov();
return true;
} else if(idname == "RNSID") {
......@@ -22,6 +22,6 @@ bool prepare_id(std::string idname) {
} else if(idname == "TOVSolver") {
id = new TovSolverAdapter();
return true;
}
}*/
return false; // no success
}
......@@ -25,9 +25,9 @@ extern idobj *id; // a global accessible pointer
bool prepare_id(std::string idname);
// C ID:
#include "InitialData/PizzaTOV.h"
#include "InitialData/RNSID.h"
#include "InitialData/TovSolverAdapter.h"
//#include "InitialData/PizzaTOV.h"
//#include "InitialData/RNSID.h"
//#include "InitialData/TovSolverAdapter.h"
// FORTRAN ID:
extern "C" {
......
......@@ -19,8 +19,8 @@ RECURSIVE SUBROUTINE InitialData(x, t, Q)
! Call InitialOrsagTang(x, 0.0 , Q)
! CALL InitialAccretionDisc(x, 0.0, Q)
!CALL InitialAccretionDisc3D(x, 0.0, Q)
CALL InitialDataTN(x, t, Q)
CALL InitialAccretionDisc3D(x, 0.0, Q)
!CALL InitialDataTN(x, t, Q)
END SUBROUTINE InitialData
......
#include "tarch/logging/Log.h"
#include "AbstractGRMHDSolver_ADERDG.h"
#include "GRMHDSolver_ADERDG_Variables.h"
......@@ -7,14 +8,14 @@
#include <stdio.h>
#ifndef PIZZATOV_AVAILABLE
/*
pizzatov::pizzatov() {
printf("Cannot call Pizza as not compiled with -DPIZZATOV_AVAILABLE");
abort();
}
void pizzatov::Interpolate(const double* x, double t, double* Q) {}
#else /* PIZZATOV_AVAILABLE */
#else
#include "pizza_tovfront/pointwise_tov.h"
......@@ -102,5 +103,5 @@ void pizzatov::Interpolate(const double* x, double t, double* Q) {
pdeprim2cons_(Q, V);
}
*/
#endif /* PIZZATOV_AVAILABLE */
/*
#include "AbstractGRMHDSolver_ADERDG.h"
#include "GRMHDSolver_ADERDG_Variables.h"
#include "InitialData/InitialData.h"
#include "Fortran/PDE.h"
#undef RNSID_AVAILABLE
#ifndef RNSID_AVAILABLE
#include <stdlib.h>
......@@ -15,7 +18,7 @@ rnsid::rnsid() {
void rnsid::Interpolate(const double* x, double t, double* Q) {}
#else /* RNSID_AVAILABLE */
#else
#include "rnsid/rnsid.h"
......@@ -112,7 +115,7 @@ void rnsid::Interpolate(const double* pos, double t, double* Q) {
} else {
get_conserved_quantities(this, pos, Q);
}
*/
/*
double V[nVar] = {0.0};
#if DIMENSIONS == 2
......@@ -130,7 +133,7 @@ void rnsid::Interpolate(const double* pos, double t, double* Q) {
#endif
*/
/*
} */
}
#endif /* RNSID_AVAILABLE */
//#endif /* RNSID_AVAILABLE */
/*
#include "TovSolverAdapter.h"
#include "Fortran/PDE.h" // P2C
......@@ -66,4 +67,4 @@ void TovSolverAdapter::Interpolate(const double* const x, double t, double* cons
pdeprim2cons_(Q, V);
}