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

Commit 70dffd3f authored by Francesco's avatar Francesco

updates GMRHDb

parent f449ee1e
......@@ -111,47 +111,46 @@ RECURSIVE SUBROUTINE PDEFluxGRMHD(F,Q)
!
F = 0.0
!
!PRINT *,"SIZE(F)",SIZE(F)
!IF( ANY( ISNAN(Q) )) THEN
! PRINT *,' Q NAN ',nVar
! PRINT *, Q
! STOP
! ENDIF
! !
!IF(ANY(Q(6:8).NE.0.0)) THEN
! PRINT *, Q(6:8)
! PRINT *, "I feel magnetized :-( PDEFluxGRMHD in"
! ERROR STOP
! ENDIF
! PRINT *, Q
! STOP
! ENDIF
! !
!IF(ANY(Q(6:8).NE.0.0)) THEN
! PRINT *, Q(6:8)
! PRINT *, "I feel magnetized :-( PDEFluxGRMHD in"
! ERROR STOP
! ENDIF
!
CALL PDECons2PrimGRMHD(V,Q,iErr)
!
CALL PDECons2PrimGRMHD(V,Q,iErr)
!
!F(1,1) = V(1)*V(2)
!RETURN
!IF(ANY(Q(6:8).NE.0.0)) THEN
! PRINT *, Q(6:8)
! PRINT *, "I feel magnetized :-( PDEFluxGRMHD in 2"
! ERROR STOP
! ENDIF
! PRINT *, Q(6:8)
! PRINT *, "I feel magnetized :-( PDEFluxGRMHD in 2"
! ERROR STOP
! ENDIF
!
gamma1 = EQN%gamma/(EQN%gamma-1.0)
rho = V(1)
!
gamma1 = EQN%gamma/(EQN%gamma-1.0)
rho = V(1)
!
!IF( ANY( ISNAN(Q) )) THEN
! PRINT *,' Q NAN ',iErr,nVar
! PRINT *, Q
! ERROR STOP
! ENDIF
! !
!IF( ANY( ISNAN(V) )) THEN
! PRINT *,' V NAN ',iErr ,nVar
! PRINT *, V
! PRINT *,' Q:'
! PRINT *, Q
! !i=20
! !V(i) = Q(1)/(Q(2)-Q(2))
! ERROR STOP
! ENDIF
! PRINT *, Q
! ERROR STOP
! ENDIF
! !
!IF( ANY( ISNAN(V) )) THEN
! PRINT *,' V NAN ',iErr ,nVar
! PRINT *, V
! PRINT *,' Q:'
! PRINT *, Q
! !i=20
! !V(i) = Q(1)/(Q(2)-Q(2))
! ERROR STOP
! ENDIF
DO i=1,3
v_cov(i) = V(1+i)
......@@ -334,7 +333,7 @@ RECURSIVE SUBROUTINE PDEFluxGRMHD(F,Q)
F(9,2) = lapse*F(9,2) - shift(2)*Q(9)
#ifdef Dim3
F(9,3) = lapse*F(9,3) - shift(3)*Q(9)
#endif
#endif
#endif
!!
!F(2:nVar,1:3) = 0.
......@@ -349,8 +348,8 @@ RECURSIVE SUBROUTINE PDEFluxGRMHD(F,Q)
! PRINT *,"uem", uem
! PRINT *,"p", p
! PRINT *, "I feel strange :-( PDEFluxGRMHD "
! ERROR STOP
! ENDIF
! ERROR STOP
! ENDIF
!
CONTINUE
!
......@@ -389,7 +388,6 @@ RECURSIVE SUBROUTINE PDEFluxPrimGRMHD(F,V,Q)
!
F = 0.0
!
PRINT *,"FLUX: nVar = ", nVar
gamma1 = EQN%gamma/(EQN%gamma-1.0)
rho = V(1)
!
......@@ -571,10 +569,10 @@ RECURSIVE SUBROUTINE PDEFluxPrimGRMHD(F,V,Q)
!
#ifdef COVCLEAN
F(9,1) = lapse*F(9,1) - shift(1)*Q(9)
F(9,2) = lapse*F(9,1) - shift(2)*Q(9)
F(9,2) = lapse*F(9,2) - shift(2)*Q(9)
#ifdef Dim3
F(9,3) = lapse*F(9,1) - shift(3)*Q(9)
#endif
F(9,3) = lapse*F(9,3) - shift(3)*Q(9)
#endif
#endif
!
CONTINUE
......@@ -819,10 +817,10 @@ RECURSIVE SUBROUTINE PDEFluxPrimVectorGRMHD(F,V,Q)
!
#ifdef COVCLEAN
F(:,9,1) = lapse(:)*F(:,9,1) - shift(:,1)*Q(:,9)
F(:,9,2) = lapse(:)*F(:,9,1) - shift(:,2)*Q(:,9)
F(:,9,2) = lapse(:)*F(:,9,2) - shift(:,2)*Q(:,9)
#ifdef Dim3
F(:,9,3) = lapse(:)*F(:,9,1) - shift(:,3)*Q(:,9)
#endif
F(:,9,3) = lapse(:)*F(:,9,3) - shift(:,3)*Q(:,9)
#endif
#endif
!
CONTINUE
......@@ -868,7 +866,7 @@ RECURSIVE SUBROUTINE PDENCPGRMHD(BgradQ,Q,gradQIn)
! ERROR STOP
!ENDIF
!
!
!
!BgradQ = 0.0
DO i=1,nVar
BgradQ(i) = 0.
......@@ -891,7 +889,7 @@ RECURSIVE SUBROUTINE PDENCPGRMHD(BgradQ,Q,gradQIn)
shift(3) = Q(13)
!
DO i=1,6
gammaij(i) = Q(13+i)
gammaij(i) = Q(13+i)
ENDDO
g_cov(1,1) = Q(14)
g_cov(1,2) = Q(15)
......@@ -1107,7 +1105,7 @@ RECURSIVE SUBROUTINE PDENCPGRMHD(BgradQ,Q,gradQIn)
!
DO i=1,nVar
BgradQ(i) = AQx(i) + BQy(i) + CQz(i)
ENDDO
ENDDO
!
IF(ABS(BgradQ(2)).gt.1e-8) THEN
......@@ -1148,7 +1146,7 @@ RECURSIVE SUBROUTINE PDENCPPrimGRMHD(BgradQ,Vc,Q,gradQ)
REAL :: AQx(nVar), BQy(nVar), CQz(nVar)
REAL :: lapse, shift(3), gammaij(6), delta(3,3), B_cov(3), vxB_cov(3), vxb_contr(3), psi, S_contr(3), qb_contr(3), B_contr(3)
REAL :: v2,v_contr(3),uem,b2,e2,gp,gm,lf,w,ww,gamma1,rho,v_cov(3), w_ij, wim
INTEGER :: i,j,k,l,m,iErr, ccount
INTEGER :: i,j,k,l,m,iErr, ccount
!
! DO dim=1,nDim
! DO i =1,nVar
......@@ -2382,7 +2380,17 @@ RECURSIVE SUBROUTINE PDEEigenvaluesGRMHD(L,Q,normal)
g_cov(2,1) = V(15)
g_cov(3,1) = V(16)
g_cov(3,2) = V(18)
!
!!
!IF(ABS(V(2)-0.2).LT.1e-4) THEN
! WRITE(*,'(a,f18.10,f18.10)') "n(1),n(2):",n(1),n(2)
! DO i = 1, nVar
! WRITE(*,'(a,i9,E16.6,E16.6,E16.6)') "i,Q(i),V(i),L(i):",i,Q(i),V(i),L(i)
! ENDDO
! WRITE(*,'(a,E16.6,E16.6,E16.6,E16.6,E16.6)') ",u,cs2,l2m,gg,den: ",u,cs2,lf2m,gg,den
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",v_cov(1),v_cov(2),v_cov(3): ",v_cov(1),v_cov(2),v_cov(3)
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",v_contr(1),v_contr(2),v_contr(3): ",v_contr(1),v_contr(2),v_contr(3)
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",n(1),n(2),n(3): ",n(1),n(2),n(3)
!ENDIF
CALL MatrixInverse3x3(g_cov,g_contr,gp)
gp = SQRT(gp)
gm = 1./gp
......@@ -2440,15 +2448,27 @@ RECURSIVE SUBROUTINE PDEEigenvaluesGRMHD(L,Q,normal)
L(i) = lapse*L(i) - sft
ENDDO
!
IF(Q(1).LT.1e-9) THEN
L = 1.0
ENDIF
!IF(Q(1).LT.1e-9) THEN
! L = 1.0
!ENDIF
!
DO i=6,nVar
L(i) = 0.
L(i) = EQN%DivCleaning_a
ENDDO
!
!
!IF(ABS(V(2)-0.2).LT.1e-4) THEN
! WRITE(*,'(a,f18.10,f18.10)') "n(1),n(2):",n(1),n(2)
! DO i = 1, nVar
! WRITE(*,'(a,i9,E16.6,E16.6,E16.6)') "i,Q(i),Vp(i),L(i):",i,Q(i),V(i),L(i)
! ENDDO
! WRITE(*,'(a,E16.6,E16.6,E16.6,E16.6,E16.6)') ",u,cs2,l2m,gg,den: ",u,cs2,lf2m,gg,den
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",v_cov(1),v_cov(2),v_cov(3): ",v_cov(1),v_cov(2),v_cov(3)
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",v_contr(1),v_contr(2),v_contr(3): ",v_contr(1),v_contr(2),v_contr(3)
! WRITE(*,'(a,E16.6,E16.6,E16.6)') ",n(1),n(2),n(3): ",n(1),n(2),n(3)
!ENDIF
!STOP
RETURN
!
! SAFE MODE: define also 'covariant' eigenvalues! (we may use the remaining free slots in L, L(6:8), L(10:19)
......@@ -2554,12 +2574,12 @@ RECURSIVE SUBROUTINE PDEEigenvectorsGRMHD(R,L,iR,Q,normal)
CALL PDECons2PrimGRMHD(VP,Q,iErr)
rho = VP(1)
DO i=1,3
v_cov(i) = V(1+i)
B_contr(i) = V(5+i)
shift(i) = V(10+i)
v_cov(i) = VP(1+i)
B_contr(i) = VP(5+i)
shift(i) = VP(10+i)
ENDDO
p = VP(5)
!
!
psi = VP(9)
lapse = VP(10)
!
......@@ -2572,7 +2592,7 @@ RECURSIVE SUBROUTINE PDEEigenvectorsGRMHD(R,L,iR,Q,normal)
g_cov(2,1) = VP(15)
g_cov(3,1) = VP(16)
g_cov(3,2) = VP(18)
!
!
CALL MatrixInverse3x3(g_cov,g_contr,gp)
gp = SQRT(gp)
gm = 1./gp
......@@ -2996,10 +3016,10 @@ RECURSIVE SUBROUTINE PDEMatrixBGRMHD(Bn,Q,normal)
A = 0.
B = 0.
C = 0.
!
ccount=0
!
DO i=1,3
!
ccount=0
!
DO i=1,3
! shift
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
j=1
......@@ -3083,9 +3103,9 @@ RECURSIVE SUBROUTINE PDEMatrixBGRMHD(Bn,Q,normal)
ENDIF
#endif
!
ENDIF
ENDDO
ENDDO
ENDIF
ENDDO
ENDDO
!lapse
!xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
j=1
......@@ -3142,7 +3162,7 @@ RECURSIVE SUBROUTINE PDEMatrixBGRMHD(Bn,Q,normal)
#endif
!
Bn = A*nv(1) + B*nv(2) + C*nv(3)
!
!
END SUBROUTINE PDEMatrixBGRMHD
......@@ -3177,8 +3197,8 @@ RECURSIVE SUBROUTINE PDECons2PrimGRMHD(V,Q,iErr)
LOGICAL :: FAILED
!
iErr = 0
tol = 1.0e-18
!
tol = 1.0e-18
!
!V = Q
!RETURN
V=0.
......
......@@ -366,25 +366,25 @@ void GRMHDb::GRMHDbSolver_ADERDG::eigenvalues(const double* const Q,const int d,
//}
//printf("\n******* ADERDG::eigenvalues *****************");
// @todo Please implement/augment if required
lambda[0] = 1.0;
lambda[1] = 1.0;
lambda[2] = 1.0;
lambda[3] = 1.0;
lambda[4] = 1.0;
lambda[5] = 1.0;
lambda[6] = 1.0;
lambda[7] = 1.0;
lambda[8] = 1.0;
lambda[9] = 1.0;
lambda[10] = 1.0;
lambda[11] = 1.0;
lambda[12] = 1.0;
lambda[13] = 1.0;
lambda[14] = 1.0;
lambda[15] = 1.0;
lambda[16] = 1.0;
lambda[17] = 1.0;
lambda[18] = 1.0;
lambda[0] = 0.0;
lambda[1] = 0.0;
lambda[2] = 0.0;
lambda[3] = 0.0;
lambda[4] = 0.0;
lambda[5] = 0.0;
lambda[6] = 0.0;
lambda[7] = 0.0;
lambda[8] = 0.0;
lambda[9] = 0.0;
lambda[10] = 0.0;
lambda[11] = 0.0;
lambda[12] = 0.0;
lambda[13] = 0.0;
lambda[14] = 0.0;
lambda[15] = 0.0;
lambda[16] = 0.0;
lambda[17] = 0.0;
lambda[18] = 0.0;
//return;
......@@ -635,14 +635,14 @@ bool GRMHDb::GRMHDbSolver_ADERDG::isPhysicallyAdmissible(
#ifdef Dim2
double dr;
dr = dx[0] * dx[0] + dx[1] * dx[1];
dr = sqrt(dr);
dr = sqrt(dr);
double radiusC;
radiusC = center[0] * center[0] + center[1] * center[1];
if (radiusC > 0.) {
radiusC = sqrt(radiusC);
}
if (radiusC + 0.5*dr > 8.05 && radiusC-0.5*dr < 8.35) {
//if (radiusC > 7.6 && radiusC < 8.5) {
radiusC = sqrt(radiusC);
}
if (radiusC + 0.5*dr > 8.05 && radiusC - 0.5*dr < 8.35) {
//if (radiusC > 7.6 && radiusC < 8.5) {
return false;
}
else {
......
......@@ -138,27 +138,35 @@ void GRMHDb::GRMHDbSolver_FV::eigenvalues(const double* const Q, const int dInde
// Dimensions = 3
// Number of variables = 19 + #parameters
//printf("\n******* EIGENVALUES FV *****************");
/*
constexpr int numberOfVariables = AbstractGRMHDbSolver_FV::NumberOfVariables;
constexpr int numberOfParameters = AbstractGRMHDbSolver_FV::NumberOfParameters;
constexpr int numberOfData = numberOfVariables + numberOfParameters;
for (int i = 0; i < numberOfData; i++) {
assertion(!std::isnan(Q[i]));
}*/
// @todo Please implement/augment if required
lambda[0] = 1.0;
lambda[1] = 1.0;
lambda[2] = 1.0;
lambda[3] = 1.0;
lambda[4] = 1.0;
lambda[5] = 1.0;
lambda[6] = 1.0;
lambda[7] = 1.0;
lambda[8] = 1.0;
lambda[9] = 1.0;
lambda[10] = 1.0;
lambda[11] = 1.0;
lambda[12] = 1.0;
lambda[13] = 1.0;
lambda[14] = 1.0;
lambda[15] = 1.0;
lambda[16] = 1.0;
lambda[17] = 1.0;
lambda[18] = 1.0;
lambda[0] = 0.0;
lambda[1] = 0.0;
lambda[2] = 0.0;
lambda[3] = 0.0;
lambda[4] = 0.0;
lambda[5] = 0.0;
lambda[6] = 0.0;
lambda[7] = 0.0;
lambda[8] = 0.0;
lambda[9] = 0.0;
lambda[10] = 0.0;
lambda[11] = 0.0;
lambda[12] = 0.0;
lambda[13] = 0.0;
lambda[14] = 0.0;
lambda[15] = 0.0;
lambda[16] = 0.0;
lambda[17] = 0.0;
lambda[18] = 0.0;
//return;
......
! GRMHDb Initial Data
#define GRMHD
#define RNSTOV
RECURSIVE SUBROUTINE PDESetup(myrank)
USE, INTRINSIC :: ISO_C_BINDING
......@@ -287,6 +285,7 @@ RECURSIVE SUBROUTINE InitialField(xGP,tGP,u0)
v0(17) = 1.0
v0(19) = 1.0
ENDIF
v0(2)=0.2
!
CASE('GRMHDAlfvenWave')
rho0 = 1.
......
......@@ -11,6 +11,12 @@ RECURSIVE SUBROUTINE PDEPrim2Cons(Q,V)
INTENT(IN) :: V
INTENT(OUT) :: Q
!
#ifdef DEBUGGONE
Q = V
STOP
RETURN
#endif
!
CALL PDEPrim2ConsGRMHD(Q,V)
!
continue
......@@ -29,6 +35,13 @@ RECURSIVE SUBROUTINE PDECons2Prim(V,Q,iErr)
INTENT(OUT) :: V
INTENT(INOUT) :: iErr
!
#ifdef DEBUGGONE
V = Q
iErr = 0
STOP
RETURN
#endif
!
CALL PDECons2PrimGRMHD(V,Q,iErr)
!
continue
......@@ -51,7 +64,7 @@ RECURSIVE SUBROUTINE PDEFlux(f,g,h,Q)
!INTENT(OUT) :: FF
! Local varialbes
INTEGER :: i
REAL :: FF(nVar,d), V(nVar)
REAL :: FF(nVar,ndim), V(nVar)
!REAL :: V(nVar)
REAL, PARAMETER :: epsilon = 1e-14
!
......@@ -262,6 +275,12 @@ RECURSIVE SUBROUTINE PDESource(S,Q)
INTENT(OUT) :: S
! --------------------------------------------
!
#ifdef DEBUGGONE
S = 0.
STOP
RETURN
#endif
!
CALL PDESourceGRMHD(S,Q)
!S = 0.
!
......@@ -355,6 +374,11 @@ RECURSIVE SUBROUTINE PDEAuxVar(aux,V,x)
! Local Variables
!
IF(nAux.GT.0) THEN
#ifdef DEBUGGONE
STOP
aux = 0.
RETURN
#endif
CALL PDEAuxVarGRMHD(aux,V,x)
ENDIF
!
......@@ -455,35 +479,36 @@ RECURSIVE SUBROUTINE HLLEMFluxFV(FL,FR,QL,QR,QavL,QavR,NormalNonZero)
fR = f1R*nv(1)+g1R*nv(2)+h1R*nv(3)
fL = f1L*nv(1)+g1L*nv(2)+h1L*nv(3)
!
IF(ANY(fR(6:8).NE.0)) THEN
PRINT *,"f1R",f1R
PRINT *,"g1R",g1R
PRINT *,"h1R",h1R
PRINT *,"f1L",f1L
PRINT *,"g1L",g1L
PRINT *,"h1L",h1L
PRINT *,"dQ",dQ
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"BEFORE: fR(6:8).NE.0",fR(6:8)
STOP
ENDIF
IF(ANY(fl(6:8).NE.0)) THEN
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"BEFORE: fl(6:8).NE.0",fl(6:8)
STOP
ENDIF
!IF(ANY(fR(6:8).NE.0)) THEN
! PRINT *,"f1R",f1R
! PRINT *,"g1R",g1R
! PRINT *,"h1R",h1R
! PRINT *,"f1L",f1L
! PRINT *,"g1L",g1L
! PRINT *,"h1L",h1L
! PRINT *,"dQ",dQ
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"BEFORE: fR(6:8).NE.0",fR(6:8)
! STOP
!ENDIF
!IF(ANY(fl(6:8).NE.0)) THEN
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"BEFORE: fl(6:8).NE.0",fl(6:8)
! STOP
!ENDIF
!USE Parameters, ONLY : d,nVar,ndim
QM = 0.5*(QL+QR)
!CALL PDECons2PrimGRMHD(VM,QM,iErr)
CALL PDEEigenvalues(LL,QL,nv)
CALL PDEEigenvalues(LR,QR,nv)
IF(ANY(QM(6:8).NE.0)) THEN
PRINT *, "HLLEMFluxFV QM(6:8)",QM(6:8)
STOP
ENDIF
!IF(ANY(QM(6:8).NE.0)) THEN
! PRINT *, "HLLEMFluxFV QM(6:8)",QM(6:8)
! STOP
!ENDIF
CALL PDEEigenvalues(LM,QM,nv)
minL = MINVAL(LL)
minM = MINVAL(LM)
......@@ -517,10 +542,9 @@ ENDIF
!
IF(QR(1).LT.1e-9.OR.QL(1).LT.1e-9) THEN
deltaL = 0.
ELSE
temp = MATMUL( deltaL, iRL )
absA = absA - sR*sL/(sR-sL)*MATMUL( RL, temp ) ! HLLEM anti-diffusion
ENDIF
ENDIF
temp = MATMUL( deltaL, iRL )
absA = absA - sR*sL/(sR-sL)*MATMUL( RL, temp ) ! HLLEM anti-diffusion
!
CALL RoeMatrix(ARoe,QL,QR,nv)
!
......@@ -545,39 +569,39 @@ ENDIF
fR = fL - Dp
fL = fL + Dm
!
IF(ANY(Dp(6:8).NE.0)) THEN
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"Dp(6:8).NE.0",Dp(6:8)
STOP
ENDIF
IF(ANY(Dm(6:8).NE.0)) THEN
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"Dm(6:8).NE.0",Dm(6:8)
STOP
ENDIF
IF(ANY(fR(6:8).NE.0)) THEN
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"fR(6:8).NE.0",fR(6:8)
STOP
ENDIF
IF(ANY(fl(6:8).NE.0)) THEN
PRINT *,"QR",QR
PRINT *,"QL",QL
PRINT *,"dQ",dQ
PRINT *,"fl(6:8).NE.0",fl(6:8)
STOP
ENDIF
! REMEMBER THE FOLLOWING: we are recursively updating qh as
! q_i^{n+1} = q_i^n - FL .... i.e. F_=F_{i+1/2}_ right flux
! q_{i+1}^{n+1} = q_i^n + FR .... i.e. FR=F_{i+1/2} left flux
! see musclhancock.cpph after "// 4. Solve Riemann problems"
!
!IF(ANY(Dp(6:8).NE.0)) THEN
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"Dp(6:8).NE.0",Dp(6:8)
! STOP
!ENDIF
!IF(ANY(Dm(6:8).NE.0)) THEN
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"Dm(6:8).NE.0",Dm(6:8)
! STOP
!ENDIF
!IF(ANY(fR(6:8).NE.0)) THEN
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"fR(6:8).NE.0",fR(6:8)
! STOP
!ENDIF
!IF(ANY(fl(6:8).NE.0)) THEN
! PRINT *,"QR",QR
! PRINT *,"QL",QL
! PRINT *,"dQ",dQ
! PRINT *,"fl(6:8).NE.0",fl(6:8)
! STOP
!ENDIF
! ! REMEMBER THE FOLLOWING: we are recursively updating qh as
! ! q_i^{n+1} = q_i^n - FL .... i.e. F_=F_{i+1/2}_ right flux
! ! q_{i+1}^{n+1} = q_i^n + FR .... i.e. FR=F_{i+1/2} left flux
! ! see musclhancock.cpph after "// 4. Solve Riemann problems"
!
END SUBROUTINE HLLEMFluxFV
......@@ -593,6 +617,13 @@ IF(ANY(Q(6:8).NE.0)) THEN
PRINT *, "PDEEigenvectors Q(6:8)",Q(6:8)
STOP
ENDIF
#ifdef DEBUGGONE
R = 0.
L = 1.0
iR = 0.
STOP
RETURN
#endif
!
CALL PDEEigenVectorsGRMHD(R,L,iR,Q,nv)
!
......
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