Commit b568c5de authored by Francesco's avatar Francesco

Merge branch 'master' of gitlab.lrz.de:exahype/ExaHyPE-Engine

parents 29af1f7c 7034224d
......@@ -2,6 +2,9 @@
*_generated.cpp
cfiles.mk
buildinfo.h
Demonstrators/**/kernels
TestCases/**/kernels
ApplicationExamples/**/kernels
# Any kind of log files created by some editors
# (swap files or temporary files)
......
......@@ -6,7 +6,7 @@
#include <cstring> // memset
#include "kernels/KernelUtils.h" // matrix indexing
#include "kernels/GaussLegendreQuadrature.h"
#include "kernels/GaussLegendreBasis.h"
constexpr int nVar = GRMHD::AbstractGRMHDSolver_ADERDG::NumberOfVariables;
......@@ -37,7 +37,7 @@ void GRMHD::GRMHDSolver_ADERDG::flux(const double* const Q,double** const F) {
pdeflux_(F[0], F[1], F[2], Q);
}
void GRMHD::GRMHDSolver_ADERDG::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int normalNonZero,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut)
void GRMHD::GRMHDSolver_ADERDG::boundaryValues(const double* const x,const double t,const double dt,const int faceIndex,const int d,const double* const fluxIn,const double* const stateIn,const double* const gradStateIn,double* const fluxOut,double* const stateOut){
double Qgp[nVar];
std::memset(stateOut, 0, nVar * sizeof(double));
std::memset(fluxOut, 0, nVar * sizeof(double));
......@@ -45,8 +45,8 @@ void GRMHD::GRMHDSolver_ADERDG::boundaryValues(const double* const x,const doubl
double F[nDim][nVar];
for(int i=0; i < basisSize; i++) { // i == time
const double weight = kernels::gaussLegendreWeights[order][i];
const double xi = kernels::gaussLegendreNodes[order][i];
const double weight = kernels::legendre::weights[order][i];
const double xi = kernels::legendre::nodes[order][i];
double ti = t + xi * dt;
initialdata_(x,&ti, Qgp);
......
......@@ -30,7 +30,7 @@ class GRMHD::GRMHDSolver_ADERDG : public GRMHD::AbstractGRMHDSolver_ADERDG {
*/
static tarch::logging::Log _log;
public:
GRMHDSolver_ADERDG(const double maximumMeshSize,const int maximumMeshDepth,const int haloCells,const int regularisedFineGridLevels,const exahype::solvers::Solver::TimeStepping timeStepping,const int DMPObservables);
GRMHDSolver_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.
......
......@@ -8,7 +8,7 @@
#include "Fortran/C2P-GRMHD.h"
#include "InitialData/InitialData.h"
#include "Fortran/MassAccretionRate.h"
#include "kernels/GaussLegendreQuadrature.h"
#include "kernels/GaussLegendreBasis.h"
#include <cmath>
#include "kernels/aderdg/generic/c/sizes.cpph"
......
......@@ -10,7 +10,7 @@
#include "C2P-MHD.h"
#include "InitialDataAdapter.h"
#include "kernels/GaussLegendreQuadrature.h"
#include "kernels/GaussLegendreBasis.h"
#include <cmath>
MHD::IntegralsWriter::IntegralsWriter(MHD::MHDSolver& solver) :
......@@ -92,11 +92,11 @@ void MHD::IntegralsWriter::mapQuantities(
//const double NumberOfUnknownsPerGridPoint = nVar;
// Gauss-Legendre weights from pos argument
double wx = kernels::gaussLegendreWeights[MHD::AbstractMHDSolver::Order][pos[0]];
double wy = kernels::gaussLegendreWeights[MHD::AbstractMHDSolver::Order][pos[1]];
double wx = kernels::legendre::weights[MHD::AbstractMHDSolver::Order][pos[0]];
double wy = kernels::legendre::weights[MHD::AbstractMHDSolver::Order][pos[1]];
double wz = 1;
#ifdef Dim3
wz = kernels::gaussLegendreWeights[MHD::AbstractMHDSolver::Order][pos[2]];
wz = kernels::legendre::weights[MHD::AbstractMHDSolver::Order][pos[2]];
#endif
// volume form for integration
......
......@@ -16,7 +16,7 @@
#include <memory>
#include <cstring>
#include <stdio.h>
#include "kernels/GaussLegendreQuadrature.h"
#include "kernels/GaussLegendreBasis.h"
#include "kernels/KernelUtils.h" // matrix indexing
tarch::logging::Log MHD::MHDSolver::_log( "MHD::MHDSolver" );
......@@ -54,8 +54,8 @@ void MHD::MHDSolver::boundaryValues(const double* const x,const double t,const d
// Integrate solution in gauss points (Qgp) in time
for(int i=0; i < basisSize; i++) { // i == time
const double weight = kernels::gaussLegendreWeights[order][i];
const double xi = kernels::gaussLegendreNodes[order][i];
const double weight = kernels::legendre::weights[order][i];
const double xi = kernels::legendre::nodes[order][i];
double ti = t + xi * dt;
alfenwave_(x, Qgp, &ti);
......
......@@ -35,6 +35,8 @@ class MHD::MHDSolver : public MHD::AbstractMHDSolver {
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
);
......
......@@ -16,9 +16,12 @@ namespace NavierStokes {
/** NEW **/
template <typename GlobalObservables>
void resetGlobalObservables(GlobalObservables& obs) {
obs.gobs(0) = -1.0;
obs.gobs(1) = -1.0;
obs.gobs(2) = 0.0;
if (obs.size() != 0) {
auto *obsRaw = obs.data();
obsRaw[0] = -1.0;
obsRaw[1] = -1.0;
obsRaw[2] = 0.0;
}
}
template <
......@@ -29,8 +32,8 @@ void mergeGlobalObservables(
GlobalObservables& obs,
ReadOnlyGlobalObservables& other) {
if (obs.size() == 0) {
return;
}
return;
}
assertion2(obs.size() == other.size(),
obs.size(), other.size());
......@@ -124,53 +127,57 @@ void mapGlobalObservablesDG(
}
}
template <typename GlobalObservables>
template <typename Solver, typename GlobalObservables>
void mapGlobalObservablesFV(
GlobalObservables& globalObservables,
const double* const luh,
Solver *solver,
GlobalObservables& globalObservables,
const double* const luh,
const tarch::la::Vector<DIMENSIONS,double>& cellSize) {
// TODO(Lukas): Please implement
}
template <typename GlobalObservables, int PatchSize, int GhostLayerWidth,
int NumberOfVariables, int NumberOfParameters>
void mapGlobalObservablesFV(
GlobalObservables& globalObservables,
const double* const Q,
const tarch::la::Vector<DIMENSIONS, double>& dx,
const std::string& scenarioName,
const PDE& ns,
const AMRSettings& amrSettings) {
// Idea: Map FV-Data to DG-Representation, and compute global variables there.
// This way, we can use the same implementation of total variance for both
// solvers. Additionally, we have the same TV before and after limiting. It is
// really slow though, so it might be a good idea to approximate TV
// differently.
// TODO(Lukas) Skip mapping if not using TV, e.g. non-gradient based crit.
if (globalObservables.size() == 0) {
return;
}
assert((PatchSize - 1) % 2 == 0);
constexpr int Order = (PatchSize - 1) / 2;
constexpr int NumberOfData = NumberOfVariables + NumberOfParameters;
// TODO(Lukas) Stack allocate this?
auto QDG =
std::vector<double>(std::pow(Order + 1, DIMENSIONS) * NumberOfData);
kernels::limiter::generic::c::projectOnDGSpace<Order + 1,
NumberOfData, GhostLayerWidth>(Q, QDG.data());
NavierStokes::mapGlobalObservablesDG(globalObservables,
QDG.data(),
dx,
scenarioName,
ns,
amrSettings,
Order,
NumberOfVariables,
NumberOfParameters);
}
if (globalObservables.size() == 0) {
return;
}
constexpr unsigned int NumberOfVariables = Solver::NumberOfVariables;
constexpr unsigned int NumberOfParameters = Solver::NumberOfParameters;
constexpr auto numberOfData = NumberOfVariables + NumberOfParameters;
constexpr unsigned int PatchSize = Solver::PatchSize;
constexpr unsigned int GhostLayerWidth = Solver::GhostLayerWidth;
constexpr auto variablesPerPatch = (PatchSize+2*GhostLayerWidth)*(PatchSize+2*GhostLayerWidth)*numberOfData;
constexpr int patchBegin = GhostLayerWidth; // patchBegin cell is inside domain
constexpr int patchEnd = patchBegin+PatchSize; // patchEnd cell is outside domain
kernels::idx3 idx_slope(PatchSize+2*GhostLayerWidth,
PatchSize+2*GhostLayerWidth,
DIMENSIONS);
double slopes[variablesPerPatch*DIMENSIONS] = {0.0};
totalVariationFV(solver, luh, cellSize, slopes);
double meanReduced = -1.0;
double varReduced = -1.0;
double n = 0.0;
for (int i = patchBegin - 1; i < patchEnd; ++i) {
for (int j = patchBegin-1; j < patchEnd; ++j) {
for (int d = 0; d < DIMENSIONS; ++d) {
const auto curTv = slopes[idx_slope(i,j,d)];
std::tie(meanReduced, varReduced) = mergeVariance(meanReduced,
curTv,
varReduced,
0.0, n,
1);
n += 1;
}
}
}
auto *rawGlobalObservables = globalObservables.data();
rawGlobalObservables[0] = meanReduced;
rawGlobalObservables[1] = varReduced;
rawGlobalObservables[2] = n;
}
} // namespace NavierStokes
......
#include "totalVariation.h"
#include "kernels/GaussLegendreBasis.h"
double totalVariation(const double* Q, int order, int numberOfVariables,
int numberOfParameters,
......@@ -13,7 +14,7 @@ double totalVariation(const double* Q, int order, int numberOfVariables,
kernels::idx4(basisSize, basisSize, basisSize, numberOfData);
#endif
const auto& quadratureWeights = kernels::gaussLegendreWeights[order];
const auto& quadratureWeights = kernels::legendre::weights[order];
double tv = 0.0;
#if DIMENSIONS == 2
// x direction (independent from the y derivatives)
......@@ -24,7 +25,7 @@ double totalVariation(const double* Q, int order, int numberOfVariables,
for (int m = 0; m < numberOfVariables; m++) {
double curGrad = 0.0;
for (int n = 0; n < basisSize; n++) { // n == matmul x
const auto t = Q[idx_lQi(k, n, m)] * kernels::dudx[order][l][n];
const auto t = Q[idx_lQi(k, n, m)] * kernels::legendre::dudx[order][l][n];
curGrad += t;
}
tv += w * std::abs(curGrad);
......@@ -41,7 +42,7 @@ double totalVariation(const double* Q, int order, int numberOfVariables,
double curGrad = 0.0;
for (int n = 0; n < basisSize; n++) { // n = matmul y
const auto t = Q[idx_lQi(n, k, m)] *
kernels::dudx[order][l][n]; /* l,n: transpose */
kernels::legendre::dudx[order][l][n]; /* l,n: transpose */
curGrad += t;
}
tv += w * std::abs(curGrad);
......@@ -60,8 +61,8 @@ double totalVariation(const double* Q, int order, int numberOfVariables,
double curGrad = 0.0;
for (int n = 0; n < basisSize; n++) {
const auto t = Q[idx_lQi(j, k, n, m)] *
kernels::dudx[order][l][n];
curGrad += t;
kernels::legendre::dudx[order][l][n];
curGrad += t;
}
tv += w * std::abs(curGrad);
}
......@@ -80,8 +81,8 @@ double totalVariation(const double* Q, int order, int numberOfVariables,
double curGrad = 0.0;
for (int n = 0; n < basisSize; n++) {
const auto t = Q[idx_lQi(j, n, k, m)] *
kernels::dudx[order][l][n];
curGrad += t;
kernels::legendre::dudx[order][l][n];
curGrad += t;