Commit 67cec377 authored by Anne Reinarz's avatar Anne Reinarz

add ADM integrals writer

parent 2e68fffe
Pipeline #171197 passed with stage
This diff is collapsed.
#ifndef __EXAHYPE_CCZ4_POSTPROCESSING_SIMPLE_ADM_INTEGRALS__
#define __EXAHYPE_CCZ4_POSTPROCESSING_SIMPLE_ADM_INTEGRALS__
#include "exahype/plotters/ascii/CSVWriter.h"
#include "peano/utils/Dimensions.h"
#include "tarch/la/Vector.h"
#include <vector>
/**
* The "simple" ADM Integrals classes implement the surface integrals on a sphere, implementing
* rather simple equations (without conformal decomposition). This is described on
* https://arxiv.org/pdf/gr-qc/0701123.pdf (or in any textbook), and I call it simple compared to
* what they suggest in terms of volume integrals and quasi-excision around punctures.
*
* Following good ExaHyPE practice, the implementation here is merely independent of ExaHyPE and
* serves as a library which can be used in an ExaHyPE plotter (here: ADMIntegralsWriter).
**/
namespace ADMIntegralsSimple {
struct Sphere {
static constexpr int sphereDim = 2;
const double radius; ///< Radius (in code units, i.e. M_Sol) where to do the surface integral
const int num_points[sphereDim]; ///< Number of equidistant grid in the 2-manifold where to integrate on
int measured_points[sphereDim]; ///< How many points where measured by a grid sweep in each dimension
static constexpr int numIntegrands = 4; /// ADM Mass scalar and Angular Momentum vector
static constexpr int idxM = 0; /// ADM Mass
static constexpr int idxJ = 1; /// ADM Angular Momentum with three quantities
double I[numIntegrands]; ///< The Integrand values determined in a grid sweep
Sphere(double radius, int ntheta=50, int nphi=50); /// < Initialize the Spherical Integration
void reset(); ///< Set Quantities and measured points to zero
void plotPatch(const double* const offsetOfPatch, const double* const sizeOfPatch, double* uPatch, double timeStamp);
void Integrand(const double* const pos, const double* const dx, const double* const Q, double* I);
double coverage() const; ///< Give a number between [0-1] which indicates the fraction of points which have been integrated
std::vector<double> times; ///< Internal bookkeeping for times in previous patches
void assertTimes(double timeStamp); ///< Fails with exception if there is local time stepping
std::string toString() const; /// Short identifier of constant data hold by this Sphere() instance
}; // struct Sphere
class MultipleSpheresWriter {
exahype::plotters::ascii::CSVWriter writer;
struct Line {
int plotindex;
double time;
double radius;
double M;
double J0, J1, J2;
double coverage;
};
Line cur_line;
public:
std::vector<Sphere> spheres;
MultipleSpheresWriter(
const std::string filenamebase,
const std::string reductionsSuffix=".asc");
// Syntactic sugar
void addSphere(double radius, int ntheta=50, int nphi=50) {
spheres.push_back(Sphere(radius,ntheta,nphi));
}
void startPlotting(double current_time);
void plotPatch(const double* const offsetOfPatch, const double* const sizeOfPatch, double* uPatch, double timeStamp);
void finishPlotting();
}; // class MultipleSpheresWriter
} // ns ADMIntegralsSimple
#endif /* __EXAHYPE_CCZ4_POSTPROCESSING_SIMPLE_ADM_INTEGRALS__ */
#include "ADMIntegralsWriter.h"
FOCCZ4::ADMIntegralsWriter::ADMIntegralsWriter() : exahype::plotters::ADERDG2UserDefined::ADERDG2UserDefined(),
simple("output/adm-integrals")
{
simple.addSphere(5);
}
void FOCCZ4::ADMIntegralsWriter::startPlotting( double time) {
simple.startPlotting(time);
}
void FOCCZ4::ADMIntegralsWriter::plotPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* u,
double timeStamp) {
simple.plotPatch(offsetOfPatch.data(), sizeOfPatch.data(), u, timeStamp);
}
void FOCCZ4::ADMIntegralsWriter::finishPlotting() {
simple.finishPlotting();
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef ADMIntegralsWriter_CLASS_HEADER_
#define ADMIntegralsWriter_CLASS_HEADER_
#include "exahype/plotters/ADERDG2UserDefined.h"
#include "ADMIntegralsSimple.h"
namespace FOCCZ4 {
class ADMIntegralsWriter;
}
class FOCCZ4::ADMIntegralsWriter : public exahype::plotters::ADERDG2UserDefined {
public:
// for the simple integrals:
ADMIntegralsSimple::MultipleSpheresWriter simple;
/**
* Constructor.
*
* \note ExaHyPE does not increment file counters for
* you if you use user defined plotting. You have
* to declare and manage such member variables yourself.
*/
ADMIntegralsWriter();
/**
* This method is invoked every time a cell
* is touched by the plotting device.
*
* \note Use the protected variables _order, _variables to
* determine the size of u.
* The array u has the size _variables * (_order+1)^DIMENSIONS.
*
* \param[in] offsetOfPatch the offset of the cell/patch.
* \param[in] sizeOfPatch the offset of the cell/patch.
* \param[in] u the degrees of freedom "living" inside of the patch.
*/
void plotPatch(
const tarch::la::Vector<DIMENSIONS, double>& offsetOfPatch,
const tarch::la::Vector<DIMENSIONS, double>& sizeOfPatch, double* u,
double timeStamp) override;
/**
* This method is called at the beginning of the plotting.
* You can use it to reset member variables, e.g., those
* used for calculations, or to increment file counters.
*
* \param[in] time a characteristic solver time stamp.
* Usually the global minimum.
*/
void startPlotting( double time) override;
/**
* This method is called at the end of the plotting.
* You can use it to reset member variables, finalise calculations (compute square roots etc.),
* or to increment file counters
*/
void finishPlotting() override;
};
#endif /* ADMIntegralsWriter_CLASS_HEADER_ */
......@@ -296,8 +296,8 @@ bool FOCCZ4::FOCCZ4Solver_ADERDG::isPhysicallyAdmissible(
{
int limvalue;
pdelimitervalue_(&limvalue,&cellCentre[0],&NumberOfDMPObservables, localDMPObservablesMin, localDMPObservablesMax);
const int num = 3;
pdelimitervalue_(&limvalue,&cellCentre[0],&num, localDMPObservablesMin, localDMPObservablesMax);
bool ret_value;
limvalue > 0 ? ret_value=false : ret_value=true;
return ret_value;
......
......@@ -304,7 +304,7 @@ SUBROUTINE NSTOV_solution(r,dr,q,printer)
dr(0) = r(1)-r(0)
q(1:NSTOV_nODE,0) = q(1:NSTOV_nODE,1)
!
IF(printer.AND.myrank.eq.0) THEN
IF(printer.eq.1.AND.myrank.eq.0) THEN
OPEN (1,file='NSTOV_q.dat')
OPEN (3,file='NSTOV_dr.dat')
100 format (1000E18.10)
......
tensish.cpph ist ein Softlink <-[mit]-> GRMHD's SVEC tensish
Grund: kann dort dann einfacher gepusht werden
This diff is collapsed.
#!/bin/bash
#
# This is a small script to copy the SVEC source code to the local application.
# We actually manage the SVEC code at https://bitbucket.org/svek/svec
# and want to keep the code there. ExaHyPE lacks a proper way of managing external
# repos...
# SvenK, 2018-07-17
SVEC="../../../ExternalLibraries/SVEC"
[[ -e $SVEC ]] || { echo "Cannot find $SVEC"; exit -1; }
svec="$SVEC/svec"
[[ -e $svec ]] || $SVEC/download-SVEC.sh || { echo "Cannot download SVEC"; exit -1; }
for remotef in $svec/src/PDE/tensish.cpph; do
localf=$(basename $remotef)
if [[ -e $localf ]] && ! diff $localf $remotef ; then
echo "$localf has local changes, compared to $remotef , please merge back!"
exit -1
else
cp -v $remotef $localf
fi
done
echo "Done"
......@@ -92,26 +92,116 @@
"implementation": "generic"
},
"variables": [
{
"name": "G",
"multiplicity": 6
},
{
"name": "MYQ",
"multiplicity": 96
"name": "K",
"multiplicity": 6
},
{
"name": "theta",
"multiplicity": 1
},
{
"name": "Z",
"multiplicity": 3
},
{
"name": "lapse",
"multiplicity": 1
},
{
"name": "shift",
"multiplicity": 3
},
{
"name": "b",
"multiplicity": 3
},
{
"name": "dLapse",
"multiplicity": 3
},
{
"name": "dxShift",
"multiplicity": 3
},
{
"name": "dyShift",
"multiplicity": 3
},
{
"name": "dzShift",
"multiplicity": 3
},
{
"name": "dxG",
"multiplicity": 6
},
{
"name": "dyG",
"multiplicity": 6
},
{
"name": "dzG",
"multiplicity": 6
},
{
"name": "traceK",
"multiplicity": 1
},
{
"name": "phi",
"multiplicity": 1
},
{
"name": "P",
"multiplicity": 3
},
{
"name": "K0",
"multiplicity": 1
},
{
"name": "rho",
"multiplicity": 1
},
{
"name": "vel",
"multiplicity": 3
},
{
"name": "pressure",
"multiplicity": 1
},
{
"name": "glm",
"multiplicity": 32
}
],
"parameters": {
"reference": "CCZ4TOV"
},
"plotters": [
{
"type": "user::defined",
"name": "TecplotWriter",
"time": 0.0,
"repeat": 1.0,
"output": "./output/tecplot",
"variables": 96
}
]
{
"type": "user::defined",
"name": "TecplotWriter",
"time": 0.0,
"repeat": 1.0,
"output": "./output/tecplot",
"variables": 96
},
{
"type": "user::defined",
"name": "ADMIntegralsWriter",
"variables": 0,
"time" : 0.0,
"repeat": 0.01,
"output": "./output/adm-integrals"
}
]
}
]
}
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