Commit 627024ae authored by Anne Reinarz's avatar Anne Reinarz

add central density, adm constraints (need debugging)

parent 67cec377
! ADM CONSTRAINTS
RECURSIVE SUBROUTINE ADMConstraints( Constraints, Q, gradQ )
USE mainVariables, ONLY : nVar, nDim, d
IMPLICIT NONE
INTENT(IN) :: Q, gradQ
INTENT(OUT) :: Constraints
INTEGER, PARAMETER :: nConstraints = 6
INTEGER :: i, ip, j, k, l, m, n, iErr, qq, ii, jj, kk, ll, mm, nn
REAL :: xGP(d), Constraints(nConstraints), Q(nVar), gradQ(nVar,nDim)!, gradQT(d,Nvar)
REAL :: traceK, R, phi, KK2
REAL :: g_contr(3,3), g_cov(3,3), Ricci(3,3)
REAL :: DD(3,3,3), Atilde(3,3), PP(3), GG(3), dP(3,3)
REAL :: s1,s2,s3,s4,s5,s6,s7,s8,s9,s10
REAL :: dDD(3,3,3,3), Christoffel(3,3,3), ChristoffelNC(3,3,3), Id(3,3), dgup(3,3,3), Riemann(3,3,3,3), dChristoffel(3,3,3,3)
REAL :: Christoffel_diff(3,3,3,3), Ham, Mom(3), dK(3,3,3), dAtilde(3,3,3), dtraceK(3), Qx(nVar), Qy(nVar), Qz(nVar)
REAL :: dg_cov(3,3,3), g_covx(3,3), g_covy(3,3), g_covz(3,3), det
REAL :: alpha, Aex(3,3), Kex(3,3), traceA, k0, dAex(3,3,3), dKex(3,3,3), Amix(3,3), Aup(3,3), Kmix(3,3), Kup(3,3)
REAL :: ghat(3), theta, dtheta(3), dghat(3,3), AA(3), dAA(3,3), BB(3,3), dBB(3,3,3), dphi(3), dPP(3,3), beta(3)
REAL :: Christoffel_tilde(3,3,3), Gtilde(3), Christoffel_kind1(3,3,3), Z(3), Zup(3), Kupdown
Constraints(:) = 0.0
Qx = gradQ(:,1)
Qy = gradQ(:,2)
if(nDim==3) then
Qz = gradQ(:,3)
else
Qz = 0
end if
!gradQT = TRANSPOSE(gradQ)
g_cov(1,1) = Q(1)
g_cov(1,2) = Q(2)
g_cov(1,3) = Q(3)
g_cov(2,1) = Q(2)
g_cov(2,2) = Q(4)
g_cov(2,3) = Q(5)
g_cov(3,1) = Q(3)
g_cov(3,2) = Q(5)
g_cov(3,3) = Q(6)
! This determinant should be close to unity, since we use the conformal decomposition
det = (Q(1)*Q(4)*Q(6)-Q(1)*Q(5)**2-Q(2)**2*Q(6)+2*Q(2)*Q(3)*Q(5)-Q(3)**2*Q(4))
g_contr(1,1) = (Q(4)*Q(6)-Q(5)**2) / det
g_contr(1,2) = -(Q(2)*Q(6)-Q(3)*Q(5))/ det
g_contr(1,3) = (Q(2)*Q(5)-Q(3)*Q(4)) / det
g_contr(2,1) = -(Q(2)*Q(6)-Q(3)*Q(5))/ det
g_contr(2,2) = (Q(1)*Q(6)-Q(3)**2) / det
g_contr(2,3) = -(Q(1)*Q(5)-Q(2)*Q(3))/ det
g_contr(3,1) = (Q(2)*Q(5)-Q(3)*Q(4)) / det
g_contr(3,2) = -(Q(1)*Q(5)-Q(2)*Q(3))/ det
g_contr(3,3) = (Q(1)*Q(4)-Q(2)**2) / det
alpha = EXP(Q(17))
!
K0 = Q(59)
!
Aex(1,1) = Q(7)
Aex(1,2) = Q(8)
Aex(1,3) = Q(9)
Aex(2,1) = Q(8)
Aex(2,2) = Q(10)
Aex(2,3) = Q(11)
Aex(3,1) = Q(9)
Aex(3,2) = Q(11)
Aex(3,3) = Q(12)
!
traceA = SUM(g_contr*Aex)
Aex = Aex - 1./3.*g_cov*traceA
!
dAex(:,1,1) = gradQ(7,:)
dAex(:,1,2) = gradQ(8,:)
dAex(:,1,3) = gradQ(9,:)
dAex(:,2,1) = gradQ(8,:)
dAex(:,2,2) = gradQ(10,:)
dAex(:,2,3) = gradQ(11,:)
dAex(:,3,1) = gradQ(9,:)
dAex(:,3,2) = gradQ(11,:)
dAex(:,3,3) = gradQ(12,:)
!
Amix = MATMUL(g_contr, Aex)
Aup = MATMUL(g_contr, Amix)
!
Theta = Q(13)
dTheta = gradQ(13,:)
!
Ghat(1) = Q(14)
Ghat(2) = Q(15)
Ghat(3) = Q(16)
dGhat(:,1) = gradQ(14,:)
dGhat(:,2) = gradQ(15,:)
dGhat(:,3) = gradQ(16,:)
!
AA(1) = Q(24)
AA(2) = Q(25)
AA(3) = Q(26)
dAA(:,1) = gradQ(24,:)
dAA(:,2) = gradQ(25,:)
dAA(:,3) = gradQ(26,:)
!
traceK = Q(54)
dtraceK = gradQ(54,:)
!
phi = EXP(Q(55))
dphi = gradQ(55,:)
PP = Q(56:58)
dPP(:,1) = gradQ(56,:)
dPP(:,2) = gradQ(57,:)
dPP(:,3) = gradQ(58,:)
!
beta(1) = Q(18)
beta(2) = Q(19)
beta(3) = Q(20)
BB(1,1) = Q(27)
BB(2,1) = Q(28)
BB(3,1) = Q(29)
BB(1,2) = Q(30)
BB(2,2) = Q(31)
BB(3,2) = Q(32)
BB(1,3) = Q(33)
BB(2,3) = Q(34)
BB(3,3) = Q(35)
!
dBB(:,1,1) = gradQ(27,:)
dBB(:,2,1) = gradQ(28,:)
dBB(:,3,1) = gradQ(29,:)
dBB(:,1,2) = gradQ(30,:)
dBB(:,2,2) = gradQ(31,:)
dBB(:,3,2) = gradQ(32,:)
dBB(:,1,3) = gradQ(33,:)
dBB(:,2,3) = gradQ(34,:)
dBB(:,3,3) = gradQ(35,:)
!
DD(1,1,1)=Q(36)
DD(1,1,2)=Q(37)
DD(1,1,3)=Q(38)
DD(1,2,1)=Q(37)
DD(1,2,2)=Q(39)
DD(1,2,3)=Q(40)
DD(1,3,1)=Q(38)
DD(1,3,2)=Q(40)
DD(1,3,3)=Q(41)
!
DD(2,1,1)=Q(42)
DD(2,1,2)=Q(43)
DD(2,1,3)=Q(44)
DD(2,2,1)=Q(43)
DD(2,2,2)=Q(45)
DD(2,2,3)=Q(46)
DD(2,3,1)=Q(44)
DD(2,3,2)=Q(46)
DD(2,3,3)=Q(47)
!
DD(3,1,1)=Q(48)
DD(3,1,2)=Q(49)
DD(3,1,3)=Q(50)
DD(3,2,1)=Q(49)
DD(3,2,2)=Q(51)
DD(3,2,3)=Q(52)
DD(3,3,1)=Q(50)
DD(3,3,2)=Q(52)
DD(3,3,3)=Q(53)
!
dDD(:,1,1,1)=gradQ(36,:)
dDD(:,1,1,2)=gradQ(37,:)
dDD(:,1,1,3)=gradQ(38,:)
dDD(:,1,2,1)=gradQ(37,:)
dDD(:,1,2,2)=gradQ(39,:)
dDD(:,1,2,3)=gradQ(40,:)
dDD(:,1,3,1)=gradQ(38,:)
dDD(:,1,3,2)=gradQ(40,:)
dDD(:,1,3,3)=gradQ(41,:)
dDD(:,2,1,1)=gradQ(42,:)
dDD(:,2,1,2)=gradQ(43,:)
dDD(:,2,1,3)=gradQ(44,:)
dDD(:,2,2,1)=gradQ(43,:)
dDD(:,2,2,2)=gradQ(45,:)
dDD(:,2,2,3)=gradQ(46,:)
dDD(:,2,3,1)=gradQ(44,:)
dDD(:,2,3,2)=gradQ(46,:)
dDD(:,2,3,3)=gradQ(47,:)
dDD(:,3,1,1)=gradQ(48,:)
dDD(:,3,1,2)=gradQ(49,:)
dDD(:,3,1,3)=gradQ(50,:)
dDD(:,3,2,1)=gradQ(49,:)
dDD(:,3,2,2)=gradQ(51,:)
dDD(:,3,2,3)=gradQ(52,:)
dDD(:,3,3,1)=gradQ(50,:)
dDD(:,3,3,2)=gradQ(52,:)
dDD(:,3,3,3)=gradQ(53,:)
!
dgup = 0.0
DO k = 1, 3
DO m = 1, 3
DO l = 1, 3
DO n = 1, 3
DO j = 1, 3
dgup(k,m,l) = dgup(k,m,l)-g_contr(m,n)*g_contr(j,l)*2*DD(k,n,j)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
Kex = Aex/phi**2 + 1./3.*traceK*g_cov/phi**2
Kmix = MATMUL( phi**2*g_contr, Kex )
Kup = MATMUL( phi**2*g_contr, Kmix )
!
Christoffel_tilde = 0.0
Christoffel = 0.0
Gtilde = 0.0
!
DO j = 1, 3
DO i = 1, 3
DO k = 1, 3
!Christoffel_kind1(i,j,k) = DD(i,j,k)+DD(j,i,k)-DD(k,i,j) ! this definition does not work !
Christoffel_kind1(i,j,k) = DD(k,i,j)+DD(j,i,k)-DD(i,j,k) ! this definition seems to work !
DO l = 1, 3
Christoffel_tilde(i,j,k) = Christoffel_tilde(i,j,k) + g_contr(k,l)*( DD(i,j,l)+DD(j,i,l)-DD(l,i,j) )
Christoffel(i,j,k) = Christoffel(i,j,k) + g_contr(k,l)*( DD(i,j,l)+DD(j,i,l)-DD(l,i,j) ) -g_contr(k,l)*( g_cov(j,l)*PP(i)+g_cov(i,l)*PP(j)-g_cov(i,j)*PP(l) )
Gtilde(i) = Gtilde(i)+2*g_contr(i,j)*g_contr(k,l)*DD(l,j,k)
ENDDO
ENDDO
ENDDO
ENDDO
Z = 0.5*MATMUL( g_cov, Ghat - Gtilde )
Zup = MATMUL(phi**2*g_contr, Z)
!
DO i = 1, 3
DO ip = 1, 3
DO m = 1, 3
DO k = 1, 3
dChristoffel(k,i,ip,m) = 0
DO l = 1, 3
dChristoffel(k,i,ip,m) = dChristoffel(k,i,ip,m) + g_contr(m,l)*(0.5*(dDD(k,i,ip,l)+dDD(i,k,ip,l))+0.5*(dDD(k,ip,i,l)+dDD(ip,k,i,l))-0.5*(dDD(k,l,i,ip)+dDD(l,k,i,ip))) &
- g_contr(m,l)*(g_cov(ip,l)*0.5*(dPP(k,i)+dPP(i,k))+g_cov(i,l)*0.5*(dPP(k,ip)+dPP(ip,k))-g_cov(i,ip)*0.5*(dPP(k,l)+dPP(l,k))) &
+dgup(k,m,l)*(DD(i,ip,l)+DD(ip,i,l)-DD(l,i,ip)) - dgup(k,m,l)*(g_cov(ip,l)*PP(i)+g_cov(i,l)*PP(ip)-g_cov(i,ip)*PP(l)) - g_contr(m,l)*( 2*DD(k,ip,l)*PP(i)+2*DD(k,i,l)*PP(ip)-2*DD(k,i,ip)*PP(l) )
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
Riemann = 0.0
DO i = 1, 3
DO ip = 1, 3
DO m = 1, 3
DO k = 1, 3
Riemann(i,k,ip,m) = dChristoffel(k,i,ip,m)-dChristoffel(ip,i,k,m)
DO j = 1, 3
Riemann(i,k,ip,m) = Riemann(i,k,ip,m) + Christoffel(i,ip,j)*Christoffel(j,k,m) - Christoffel(i,k,j)*Christoffel(j,ip,m)
ENDDO
ENDDO
ENDDO
ENDDO
ENDDO
!
Ricci = 0.0
DO m = 1, 3
DO n = 1, 3
DO l = 1, 3
Ricci(m,n) = Ricci(m,n) + Riemann(m,l,n,l)
ENDDO
ENDDO
ENDDO
!
R = SUM(phi**2*g_contr*Ricci)
!
Kupdown = SUM(Kex*Kup)
Ham = R - KupDown + traceK**2
!
dKex = 0.0
DO j = 1, 3
DO i = 1, 3
DO k = 1, 3
dKex(k,i,j) = 1.0/phi**2*( dAex(k,i,j) - 2.0*Aex(i,j)*PP(k) + 1./3.*dtraceK(k)*g_cov(i,j) + 2./3.*traceK*DD(k,i,j) - 2./3.*traceK*g_cov(i,j)*PP(k) )
ENDDO
ENDDO
ENDDO
!
Mom(:) = 0.0
DO ii = 1, 3
DO jj = 1, 3
DO ll = 1, 3
Mom(ii) = Mom(ii) + g_contr(jj,ll)*(dKex(ll,ii,jj) - dKex(ii,jj,ll))
DO mm = 1, 3
!Mom(ii) = Mom(ii) + g_contr(jj,ll)*( - ChristoffelNC(jj,ll,mm)*K(mm,ii) + ChristoffelNC(jj,ii,mm)*K(mm,ll))
Mom(ii) = Mom(ii) + g_contr(jj,ll)*( - Christoffel(jj,ll,mm)*Kex(mm,ii) + Christoffel(jj,ii,mm)*Kex(mm,ll))
ENDDO
ENDDO
ENDDO
ENDDO
Constraints(1) = Ham
Constraints(2:4) = Mom(1:3)
Constraints(5) = det - 1.0
Constraints(6) = traceA
!
CONTINUE
!
END SUBROUTINE ADMConstraints
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
//
// ========================
// www.exahype.eu
// ========================
#include "CentralDensity.h"
FOCCZ4::CentralDensity::CentralDensity(FOCCZ4::FOCCZ4Solver& solver) {
// @TODO Please insert your code here.
}
FOCCZ4::CentralDensity::~CentralDensity() {
}
void FOCCZ4::CentralDensity::startPlotting( double time) {
// @TODO Please insert your code here.
}
void FOCCZ4::CentralDensity::finishPlotting() {
// @TODO Please insert your code here.
}
void FOCCZ4::CentralDensity::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* const Q,
double* const outputQuantities,
double timeStamp
) {
const int writtenUnknowns = 1;
outputQuantities[0] = Q[60];
}
// This file was generated by the ExaHyPE toolkit.
// It will not be overwritten.
//
// ========================
// www.exahype.eu
// ========================
#ifndef POSTPROCESSING_CentralDensity_CLASS_HEADER_
#define POSTPROCESSING_CentralDensity_CLASS_HEADER_
#include "exahype/plotters/Plotter.h"
namespace FOCCZ4 {
class FOCCZ4Solver;
class CentralDensity;
}
class FOCCZ4::CentralDensity : public exahype::plotters::Plotter::UserOnTheFlyPostProcessing {
public:
CentralDensity(FOCCZ4::FOCCZ4Solver& solver);
virtual ~CentralDensity();
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* const Q,
double* const outputQuantities,
double timeStamp) override;
};
#endif /* POSTPROCESSING_CentralDensity_CLASS_HEADER_ */
\ No newline at end of file
#include "ConstraintsWriter.h"
#include "kernels/GaussLegendreBasis.h"
#include "kernels/KernelUtils.h"
#include "PDE.h" // ADMConstraints()
#include "AbstractFOCCZ4Solver_ADERDG.h"
#include "kernels/aderdg/generic/c/computeGradients.cpph"
#include "kernels/GaussLegendreBasis.h"
#include "tarch/logging/Log.h"
static tarch::logging::Log _log("FOCCZ4::ConstraintsWriter::");
FOCCZ4::ConstraintsWriter::ConstraintsWriter(exahype::solvers::LimitingADERDGSolver& solver) {
// @todo Please insert your code here
}
FOCCZ4::ConstraintsWriter::ConstraintsWriter(FOCCZ4Solver_FV& solver) {
// @todo Please insert your code here
}
FOCCZ4::ConstraintsWriter::ConstraintsWriter(FOCCZ4Solver_ADERDG& solver) {
// @todo Please insert your code here
}
FOCCZ4::ConstraintsWriter::~ConstraintsWriter() {
// @todo Please insert your code here
}
void FOCCZ4::ConstraintsWriter::startPlotting(double time) {
// @todo Please insert your code here
timeLastWarned = -1;
}
void FOCCZ4::ConstraintsWriter::finishPlotting() {
// @todo Please insert your code here
}
bool FOCCZ4::ConstraintsWriter::mapWithDerivatives() {
return true;
}
void FOCCZ4::ConstraintsWriter::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* gradQ,
double* outputQuantities,
double timeStamp
) {
// this should be a UserOnTheFlyPostProcessing constant,
// allowing to ensure we write out 6 unknowns.
static constexpr int writtenUnknowns = 6;
for(int i=0; i<59; i++) { if(Q[i]!=Q[i]) std::abort(); }
admconstraints_(outputQuantities, Q, gradQ);
for(int i=0; i<6; i++) { if(outputQuantities[i]!=outputQuantities[i]) std::abort(); }
}
void FOCCZ4::ConstraintsWriter::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) {
// this should be a UserOnTheFlyPostProcessing constant,
// allowing to ensure we write out 6 unknowns.
static constexpr int writtenUnknowns = 6;
//for(int i=0; i<59; i++) { if(Q[i]!=Q[i]) std::abort(); }
double gradQ[basisSize3 * DIMENSIONS * numberOfVariables];
//kernels::aderdg::generic::c::computeGradQ<FOCCZ4::AbstractFOCCZ4Solver_ADERDG>(gradQ, Q, sizeOfPatch);
admconstraints_(outputQuantities, Q, gradQ);
for(int i=0; i<6; i++) { if(outputQuantities[i]!=outputQuantities[i]) std::abort(); }
}
// 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 "exahype/plotters/Plotter.h"
#include "exahype/solvers/LimitingADERDGSolver.h"
#include "AbstractFOCCZ4Solver_ADERDG.h"
namespace FOCCZ4{
class ConstraintsWriter;
/**
* Forward declaration
*/
class FOCCZ4Solver_ADERDG;
class FOCCZ4Solver_FV;
}
class FOCCZ4::ConstraintsWriter: public exahype::plotters::Plotter::UserOnTheFlyPostProcessing{
double timeLastWarned;
public:
static constexpr int numberOfVariables = FOCCZ4::AbstractFOCCZ4Solver_ADERDG::NumberOfVariables;
static constexpr int order = FOCCZ4::AbstractFOCCZ4Solver_ADERDG::Order;
static constexpr int basisSize = order + 1;
static constexpr int basisSize3 = basisSize*basisSize*basisSize;
public:
ConstraintsWriter(FOCCZ4Solver_ADERDG& solver);
ConstraintsWriter(FOCCZ4Solver_FV& solver);
ConstraintsWriter(exahype::solvers::LimitingADERDGSolver& solver);
virtual ~ConstraintsWriter();
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* gradQ,
double* outputQuantities,
double timeStamp) 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;
bool mapWithDerivatives() override;
};
......@@ -29,7 +29,9 @@ void pderefinecriteria_(int* refine_flag,const double* max_luh,const double* min
//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);
void admconstraints_(double* constraints, double* Q, double* gradQ);
}/* extern "C" */
#endif /* __EXAHYPE_USER_PDE__ */
......@@ -88,7 +88,7 @@
"source"
],
"scheme": "musclhancock",
"slope_limiter" : "minmod",
"slope_limiter" : "mclim",
"implementation": "generic"
},
"variables": [
......@@ -200,6 +200,27 @@
"time" : 0.0,
"repeat": 0.01,
"output": "./output/adm-integrals"
},
{
"type": "vtk::Cartesian::vertices::binary",
"name": "ConstraintsWriter",
"variables": 6,
"time": 0.0,
"repeat": 0.01,
"output": "./vtk-output/constraints"
},
{
"type": "probe::ascii",
"name": "CentralDensity",
"time": 0.0,
"repeat": 0.001,
"output": "./output/rho",
"variables": 1,
"select": {
"x": 0.0,
"y": 0.0,
"z": 0.0
}
}
]
}
......
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