Commit 7034224d authored by m.tavelli's avatar m.tavelli

First running version

parent ff171897
Pipeline #168187 failed with stage
......@@ -42,82 +42,7 @@ RECURSIVE SUBROUTINE PDEPrim2Cons(Q,V)
!
REAL :: A(3,3), G(3,3), Id(3,3), devG(3,3), eh, evv, T, falpha, Temp(3,3)
!
#ifdef EULER
!
Q(1) = V(1) ! fluid density
Q(2:4) = V(1)*V(2:4) ! momentum
Q(5) = V(5)/(EQN%gamma-1) + 0.5*V(1)*SUM(V(2:4)**2) ! total energy = internal energy + kinetic energy
!
#ifdef NUMENTR
Q(6) = V(6)
#endif
#ifdef BURGER
Q(1:6) = V(1:6)
#endif
#endif
!
#ifdef HSM
Q(1) = V(1) ! water depth
Q(2:5) = V(1)*V(2:5) ! momentum and aux. variables
Q(6) = V(6) ! still water depth
#endif
!
#ifdef HGNB
Q(1) = V(1) ! water depth
Q(2:7) = V(1)*V(2:7) ! momentum, pressure, aux. variables
Q(8) = V(8) ! bottom
#endif
!
#ifdef MAXWELL
Q(:) = V(:)
#endif
!
#ifdef ECHGNB
Q(1) = V(1) ! water depth
Q(2:7) = V(1)*V(2:7) ! momentum, sigma, p, pb
Q(8) = V(8) ! bottom
#endif
!
#ifdef HSCHROEDINGER
Q = V
Q(1) = V(1)
Q(2) = V(1)*V(2)
Q(3) = V(1)*V(3)
Q(4) = V(1)*V(4)
Q(5) = V(1)*V(5)
Q(6) = V(1)*V(6)
#endif
!
#ifdef THETAMODEL
Q(1) = V(1) ! h
Q(2) = V(1)*V(4)*V(2) ! h u
Q(3) = V(1)*V(4)*V(3) ! h v
Q(4) = V(1)*V(4) ! h theta
Q(5) = V(5) ! b
#endif
!
#ifdef ELASTICITY
Q = V
#endif
!
#ifdef ACOUSTIC
Q = V
#endif
!
#if defined(Z4EINSTEIN) || defined(Z4GRMHD)
Q = V
!
#ifdef Z4GRMHD
VGRMHD(1:9) = V(55:63) ! hydro variables
VGRMHD(10) = V(17) ! lapse
VGRMHD(11:13) = V(18:20) ! shift
VGRMHD(14:19) = V(1:6) ! metric
CALL PDEPrim2ConsGRMHD(QGRMHD,VGRMHD)
Q(55:63) = QGRMHD(1:9)
#endif
!
#endif
!
#if defined(CCZ4EINSTEIN) || defined(CCZ4GRMHD) || defined(CCZ4GRHD) || defined(CCZ4GRGPR)
Q = V
Q(17) = LOG(V(17))
......
......@@ -188,13 +188,14 @@ class FOCCZ4::FOCCZ4Solver_ADERDG : public FOCCZ4::AbstractFOCCZ4Solver_ADERDG {
void static fusedSource(const double* const restrict Q, const double* const restrict gradQ, double* const restrict S);
bool isPhysicallyAdmissible(
const double* const solution,
const double* const observablesMin,const double* const observablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& center,
const tarch::la::Vector<DIMENSIONS,double>& dx,
const double t) const override;
bool isPhysicallyAdmissible(
const double* const solution,
const double* const localDMPObservablesMin,
const double* const localDMPObservablesMax,
const bool wasTroubledInPreviousTimeStep,
const tarch::la::Vector<DIMENSIONS,double>& cellCentre,
const tarch::la::Vector<DIMENSIONS,double>& cellSize,
const double timeStamp) const override;
/* pointSource() function not included, as requested in the specification file */
/* multiplyMaterialParameterMatrix() not included, as requested in the specification file */
......
......@@ -15,10 +15,31 @@ RECURSIVE SUBROUTINE InitParameters(STRLEN,PARSETUP)
print *, 'Chosen setup: ',ICType
print *, "****************************************************************"
!stop
ICType = 'CCZ4MinkowskiSrc'
select case(ICType)
case('NLOPRUPTURE')
case('CCZ4MinkowskiSrc')
EQN%CCZ4GLMc0 = 1.5 ! 0.1
EQN%CCZ4GLMc = 1.2 ! 2.0
EQN%CCZ4GLMd = 2.0 ! 1.0
EQN%CCZ4GLMepsA = 1.0 ! 5.
EQN%CCZ4GLMepsP = 1.0 ! 5.
EQN%CCZ4GLMepsD = 1.0 ! 0.1
!
EQN%CCZ4k1 = 0.0
EQN%CCZ4k2 = 0.0
EQN%CCZ4k3 = 0.0
EQN%CCZ4eta = 0.0
EQN%CCZ4f = 0.0
EQN%CCZ4g = 0.0
EQN%CCZ4xi = 0.0
EQN%CCZ4e = 2.0
EQN%CCZ4c = 0.0
EQN%CCZ4mu = 0.0
EQN%CCZ4ds = 1.0
EQN%CCZ4sk = 0.0
EQN%CCZ4bs = 0.0 ! set bs=1 if you want to activate the shift convection for beta, b and B (standard CCZ4 formulation). set it to bs=0 to switch off shift convection for those quantities
EQN%CCZ4LapseType = 0 ! harmonic lapse
EQN%EinsteinAutoAux = 0
end select
END SUBROUTINE InitParameters
......@@ -29,17 +50,42 @@ RECURSIVE SUBROUTINE InitialData(xGP, tGp, Q)
! Argument list
REAL, INTENT(IN) :: xGP(3), tGp !
REAL, INTENT(OUT) :: Q(nVar) !
REAL :: up(nVar)
REAL :: up(nVar),V0(nVar),r
select case(ICType)
case('NLOPRUPTURE')
case('CCZ4MinkowskiSrc')
!
V0(:) = 0.0
!
V0(1) = 1.0
V0(4) = 1.0
V0(6) = 1.0
V0(17) = 1.0
V0(55) = 1.0
!
!V0(54) = 0.2*EXP(-0.5*( xGP(1)**2+xGP(2)**2)/1.0**2 )
!
r = SQRT( xGP(1)**2 + xGP(2)**2 )
!
V0(60) = 0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
IF( r>5.0 .AND. r<10. ) THEN
V0(61) = -0.2*xGP(2)*( 10.0 - r )/5.0
V0(62) = +0.2*xGP(1)*( 10.0 - r )/5.0
ELSEIF( r<=5.0 ) THEN
V0(61) = -0.2*xGP(2)
V0(62) = +0.2*xGP(1)
ELSE
V0(61:62) = 0.0
ENDIF
V0(63) = 0.0
V0(64) = 0.5*0.001*EXP(-0.5*( (xGP(1)-2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 ) + 0.5*0.001*EXP(-0.5*( (xGP(1)+2.0)**2+xGP(2)**2+0*xGP(3)**2)/1.0**2 )
!
CASE DEFAULT
PRINT *, 'NOT IMPLEMENTED'
STOP
END SELECT
CALL PDEPrim2Cons(Q,up)
CALL PDEPrim2Cons(Q,V0)
END SUBROUTINE InitialData
RECURSIVE SUBROUTINE PDElimitervalue(limiter_value,xx,numberOfObservables, observablesMin, observablesMax)
......@@ -55,18 +101,6 @@ RECURSIVE SUBROUTINE PDElimitervalue(limiter_value,xx,numberOfObservables, obser
limiter_value=0
IF((observablesMax(1)<0.999 .and. observablesMin(1)>0.001) .or. observablesMax(1)>1.001) THEN ! was -0.001 .or. levelmin<0.001
limiter_value = 1
ENDIF
IF(observablesMax(2)>1.e-5 .and. observablesMax(1)>1.e-4) THEN
limiter_value = 1
ENDIF
IF(observablesMax(3)>0.5 ) THEN
limiter_value = 1
ENDIF
!limiter_value = 1
END SUBROUTINE PDElimitervalue
......
......@@ -12,7 +12,7 @@ void pdeflux_(double* Fx, double* Fy, double* Fz, const double* const Q);
void pdesource_(double* S, const double* const Q);
void pdencp_(double* BgradQ, const double* const Q, const double* const gradQ);
void pdefusedsrcncp_(double* S, const double* const Q, const double* const gradQ)
void pdefusedsrcncp_(double* S, const double* const Q, const double* const gradQ);
void pdeeigenvalues_(double* lambda, const double* const Q, double* nv);
void pdevarname_(char* MyNameOUT, int* ind);
......
......@@ -96,11 +96,19 @@
}
],
"parameters": {
"reference": "NLOPRUPTUREPW"
"reference": "CCZ4MinkowskiSrc"
},
"plotters": [
{
"type": "user::defined",
"name": "TecplotWriter",
"time": 0.0,
"repeat": 0.1,
"output": "./output/tecplot",
"variables": 96
}
]
}
]
......
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