Dear all,
I'm trying to solve a wave equation that describes radial oscillations
of a TOV star, depending on radius r and t in Llama multipatch
coordinates (Thornburg04).
The equation has a structure of
Ẍ=AX′+BX″+CX\ddot{X} = A X' + B X'' + C X
I rewrote it as set of coupled equations to evolve it with the MoL thorn:
X˙=Π
H˙=Π′\dot{H} = \Pi'
Π˙=AH+BH′+CX\dot{\Pi} = AH + B H' +C X
or as noted in my code as it is attached:
Pidot = A*H + B*drH + C*Xi
Hdot = drPi
Xidot = Pi
For MoL registered as evolved variables are Xi with rhs Xidot, Pi with
rhs Pidot and H with rhs Hdot.
The spatial derivatives drH and drPi are also calculated during the
schedule of MoL_CalcRHS in my thorn.
Boundary conditions are applied for r>R and an interval near the TOV
radius r=R.
if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <=
(TOV_surface + rprec) )then
H(i,j,k) = 0
drPi(i,j,k) = 0
else if ( grid_r(i,j,k) > (TOV_surface + rprec) ) then
Xi(i,j,k) = 0
Pi(i,j,k) = 0
H(i,j,k) = 0
end if
if (grid_r(i,j,k) > (TOV_surface + rprec) )then
Xidot(i,j,k) = 0
Hdot(i,j,k) = 0
Pidot(i,j,k) = 0
else if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND.
grid_r(i,j,k) <= (TOV_surface + rprec) )then
Hdot(i,j,k) = 0
Pidot(i,j,k) = B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)
end if
Problems occur close to the surface of the star. My evolved variables
start do diverge close to the surface after a few iterations. Looking at
the data it seems that the divergence is founded by the values of H. I
tried following things to encircle the issues:
* If I put Pidot = B*drH + C*Xi the values seem to be ok, whereas for
Pidot = A*H the divergence appears. The A-factor only amplifies this
behavior. If I put Pidot = H it behaves the same, but much slower.
* If I use drXi instead of H (as it is commented out), it does not
make any difference.
* If I enlarge the size of the interval (e.g. TOV_surface + 5*rprec)
the same divergence appears, but shifted towards grid points next to
the interval. Also the divergence appears a few iterations later.
* If I change the drXi or H values close to the surface manually (e.g.
using a backsided differentiation, or put specific values by hand)
it also shifts the divergence (like above).
So far I don't know how to remedy this issue, but maybe I'm overlooking
something obvious.
Does anyone have an idea on what I could try?
Thanks a lot!
Best regards and merry Christmas,
Severin Frank
#===============================================================
# Testing SWTNS on a static TOV star
#
# Severin Frank
#===============================================================
#################### Thorns ####################################
ActiveThorns = "Time MoL"
ActiveThorns = "TOVSolver ADMBase HydroBase Constants StaticConformal SpaceMask"
ActiveThorns = "CoordBase CartGrid3D SymBase InitBase Boundary"
ActiveThorns = "LocalInterp AEILocalInterp LoopControl Carpet"
ActiveThorns = "IOUtil CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic"
ActiveThorns = "NaNChecker CarpetReduce"
ActiveThorns = "Coordinates Interpolate2 CarpetInterp2 GlobalDerivative SWTNS"
#CoordinatesSymmetry
#################### Terminate #################################
Cactus::terminate = "time"
Cactus::cctk_final_time = 2
#################### Time Integration ##########################
Time::timestep_method = "given"
Time::timestep = 0.09
MoL::ODE_Method = "Generic"
MoL::Generic_Type = "RK"
MoL::MoL_Intermediate_Steps = 4
MoL::MoL_Num_Scratch_Levels = 3
MoL::verbose = "extreme"
#################### Grid Setup ################################
Carpet::domain_from_multipatch = yes
CartGrid3D::type = "multipatch"
CartGrid3D::set_coordinate_ranges_on = "all maps"
Coordinates::coordinate_system = "Thornburg04"
Coordinates::h_cartesian = 0.1
Coordinates::h_radial = 0.1
SWTNS::rprec = 0.1
TOVSolver::tovrprec = 0.1
SWTNS::SWTNS_amplitude = 0.001 #0.0031830989
Coordinates::sphere_inner_radius = 2 #
Coordinates::sphere_outer_radius = 8 #7.5284600798562327 #100
Coordinates::n_angular = 9
Driver::ghost_size = 3
Coordinates::patch_boundary_size = 3
Coordinates::outer_boundary_size = 3
Coordinates::additional_overlap_size = 3
#stagger_patch_boundaries
Interpolate2::interpolator_order = 4
Interpolate2::continue_if_selftest_fails = no
#Interpolate2::verbose = yes
Coordinates::store_inverse_jacobian = "yes"
#################### Initial Data ##############################
InitBase::initial_data_setup_method = "init_some_levels"
HydroBase::timelevels = 1
ADMBase::metric_type = "physical"
ADMBase::initial_data = "tov"
ADMBase::initial_lapse = "tov"
ADMBase::initial_shift = "tov"
ADMBase::initial_dtlapse = "zero"
ADMBase::initial_dtshift = "zero"
SpaceMask::use_mask = "yes"
TOVSolver::TOV_Rho_Central[0] = 1.62e-3 #1e-3 1.62e-3 3.24e-3 1.24e-3
TOVSolver::TOV_Gamma = 2.0
SWTNS::TOV_Gamma = 2.0
TOVSolver::TOV_K = 100.0
SWTNS::TOV_K = 100.0
#TOVSolver::TOV_Num_Radial = 10000
#################### Finite Differencing ########################
SummationByParts::order = 2
SummationByParts::use_dissipation = no
SummationByParts::onesided_interpatch_boundaries = no
SummationByParts::onesided_outer_boundaries = yes
SummationByParts::sbp_upwind_deriv = yes
SummationByParts::sbp_1st_deriv = yes
SummationByParts::sbp_2nd_deriv = yes
#################### NaN Checker ################################
NaNChecker::check_every = 1
NaNChecker::action_if_found = "just warn" #"terminate", "just warn", "abort"
NaNChecker::check_vars = "HydroBase::rho HydroBase::press TOVSolver::TOV_mr
SWTNS::TOV_surface SWTNS::PHI SWTNS::drPHI
SWTNS::LAMBDA
SWTNS::drLAMBDA SWTNS::drpress SWTNS::drrho
SWTNS::A SWTNS::B SWTNS::C SWTNS::omega
SWTNS::Hdot SWTNS::drH SWTNS::drPi SWTNS::H
SWTNS::Xi SWTNS::Xidot SWTNS::Pi SWTNS::Pidot
SWTNS::Xeta #SWTNS::drXi SWTNS::param_dr
SWTNS::deriv_error"
#################### Output #####################################
IOBasic::outInfo_every = 1
IOBasic::outInfo_vars = "
#HydroBase::rho
#HydroBase::press
#SWTNS::gappa
SWTNS::Xi
#SWTNS::Xeta
#SWTNS::drXi
SWTNS::H
SWTNS::drH
#SWTNS::omega
#SWTNS::deriv_error
#SWTNS::TOV_surface
#TOVSolver::TOV_mr
#SWTNS::drPHI
#SWTNS::Xidot
#SWTNS::Pi
#SWTNS::Pidot
#SWTNS::PHI
#SWTNS::LAMBDA
#TOVSolver::TOV_PHI
#TOVSolver::TOV_R
"
#CarpetIOASCII::out1d_vars = "SWTNS::Xi"
#CarpetIOASCII::out1d_z = yes
IOASCII::compact_format = yes
# IOASCII::use_grid_coordinates = yes
#IOASCII::separate_grids = no
# IOASCII::one_file_per_group =yes
# IOASCII::out0D_every = 1
# IOASCII::out0D_point_x = 0
# IOASCII::out0D_point_y = 0
# IOASCII::out0D_point_z = 2
# IOASCII::out0D_vars ="
#
grid::coordinates{out_every=1000000000}
# SWTNS::Xeta
# SWTNS::Xi
# SWTNS::TOV_surface
#
#SWTNS::PHI{out_every=100000000000000}
#
#SWTNS::LAMBDA{out_every=100000000000000}
# SWTNS::Xidot
# SWTNS::Pi
# SWTNS::Pidot
#
#HydroBase::rho{out_every=1000000000}
#
#HydroBase::press{out_every=1000000000}
# "
IOHDF5::out_every = 1
IOHDF5::out_vars = "
grid::coordinates{out_every=100000000000}
Coordinates::patchnumber{out_every=100000000000000}
HydroBase::rho{out_every=1000000000000000}
HydroBase::press{out_every=100000000000000}
SWTNS::A{out_every=100000000000000}
SWTNS::B{out_every=100000000000000}
SWTNS::C{out_every=100000000000000}
SWTNS::omega{out_every=100000000000000}
SWTNS::drpress{out_every=100000000000000}
SWTNS::drrho{out_every=100000000000000}
SWTNS::dxdr{out_every=100000000000000}
SWTNS::dydr{out_every=100000000000000}
SWTNS::dzdr{out_every=100000000000000}
SWTNS::Xi
SWTNS::Xeta
SWTNS::Xidot
SWTNS::Pi
SWTNS::drPi
SWTNS::Pidot
SWTNS::drXi
SWTNS::drH
SWTNS::H
SWTNS::Hdot
SWTNS::grid_r{out_every=1000000000000000}
SWTNS::param_dr{out_every=1000000000000000}
SWTNS::deriv_error{out_every=1000000000000000}
TOVSolver::TOV_mr{out_every=100000000000000}
SWTNS::PHI{out_every=100000000000000}
SWTNS::drPHI{out_every=100000000000000}
SWTNS::LAMBDA{out_every=100000000000000}
SWTNS::drLAMBDA{out_every=100000000000000}
"
# This Thorn is based on LlamaWaveToy and SevesWaveToy
# Interface definition for thorn SWTNS
IMPLEMENTS: SWTNS
INHERITS: grid, Coordinates, GlobalDerivative, SummationByParts, Interpolate,
ADMBase, HydroBase, Constants, StaticConformal, TOVSolver
USES INCLUDE: constants.h
CCTK_INT FUNCTION \
MultiPatch_GetMap \
(CCTK_POINTER_TO_CONST IN cctkGH)
REQUIRES FUNCTION MultiPatch_GetMap
CCTK_INT FUNCTION \
MultiPatch_GetBbox \
(CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN size, \
CCTK_INT OUT ARRAY bbox)
REQUIRES FUNCTION MultiPatch_GetBbox
SUBROUTINE globalDiff_gv ( CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN dir, \
CCTK_REAL IN ARRAY var, \
CCTK_REAL OUT ARRAY dvar, \
CCTK_REAL IN ARRAY J_dxda, \
CCTK_REAL IN ARRAY J_dxdb, \
CCTK_REAL IN ARRAY J_dxdc, \
CCTK_REAL IN ARRAY J_dyda, \
CCTK_REAL IN ARRAY J_dydb, \
CCTK_REAL IN ARRAY J_dydc, \
CCTK_REAL IN ARRAY J_dzda, \
CCTK_REAL IN ARRAY J_dzdb, \
CCTK_REAL IN ARRAY J_dzdc, \
CCTK_INT IN table_handle )
USES FUNCTION globalDiff_gv
SUBROUTINE globalDiff2_gv( CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN dir1, \
CCTK_INT IN dir2, \
CCTK_REAL IN ARRAY var, \
CCTK_REAL OUT ARRAY dvar, \
CCTK_REAL IN ARRAY J_dadx, \
CCTK_REAL IN ARRAY J_dbdx, \
CCTK_REAL IN ARRAY J_dcdx, \
CCTK_REAL IN ARRAY J_dady, \
CCTK_REAL IN ARRAY J_dbdy, \
CCTK_REAL IN ARRAY J_dcdy, \
CCTK_REAL IN ARRAY J_dadz, \
CCTK_REAL IN ARRAY J_dbdz, \
CCTK_REAL IN ARRAY J_dcdz, \
CCTK_REAL IN ARRAY dJ_dadxdx, \
CCTK_REAL IN ARRAY dJ_dbdxdx, \
CCTK_REAL IN ARRAY dJ_dcdxdx, \
CCTK_REAL IN ARRAY dJ_dadxdy, \
CCTK_REAL IN ARRAY dJ_dbdxdy, \
CCTK_REAL IN ARRAY dJ_dcdxdy, \
CCTK_REAL IN ARRAY dJ_dadxdz, \
CCTK_REAL IN ARRAY dJ_dbdxdz, \
CCTK_REAL IN ARRAY dJ_dcdxdz, \
CCTK_REAL IN ARRAY dJ_dadydy, \
CCTK_REAL IN ARRAY dJ_dbdydy, \
CCTK_REAL IN ARRAY dJ_dcdydy, \
CCTK_REAL IN ARRAY dJ_dadydz, \
CCTK_REAL IN ARRAY dJ_dbdydz, \
CCTK_REAL IN ARRAY dJ_dcdydz, \
CCTK_REAL IN ARRAY dJ_dadzdz, \
CCTK_REAL IN ARRAY dJ_dbdzdz, \
CCTK_REAL IN ARRAY dJ_dcdzdz, \
CCTK_INT IN table_handle )
USES FUNCTION globalDiff2_gv
CCTK_INT FUNCTION \
MoLRegisterEvolvedGroup \
(CCTK_INT IN EvolvedIndex, \
CCTK_INT IN RHSIndex)
REQUIRES FUNCTION MoLRegisterEvolvedGroup
CCTK_INT FUNCTION \
MoLRegisterEvolved \
(CCTK_INT IN EvolvedIndex, \
CCTK_INT IN RHSIndex)
REQUIRES FUNCTION MoLRegisterEvolved
CCTK_INT FUNCTION \
Boundary_SelectGroupForBC \
(CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN faces, \
CCTK_INT IN boundary_width, \
CCTK_INT IN table_handle, \
CCTK_STRING IN group_name, \
CCTK_STRING IN bc_name)
REQUIRES FUNCTION Boundary_SelectGroupForBC
CCTK_INT FUNCTION \
Boundary_SelectVarForBC \
(CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN faces, \
CCTK_INT IN boundary_width, \
CCTK_INT IN table_handle, \
CCTK_STRING IN var_name, \
CCTK_STRING IN bc_name)
REQUIRES FUNCTION Boundary_SelectVarForBC
CCTK_INT FUNCTION \
Boundary_SelectedGVs \
(CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN array_size, \
CCTK_INT ARRAY OUT var_indicies, \
CCTK_INT ARRAY OUT faces, \
CCTK_INT ARRAY OUT boundary_widths, \
CCTK_INT ARRAY OUT table_handles, \
CCTK_STRING IN bc_name)
REQUIRES FUNCTION Boundary_SelectedGVs
CCTK_INT FUNCTION \
SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
REQUIRES FUNCTION SymmetryTableHandleForGrid
CCTK_POINTER_TO_CONST FUNCTION SymmetryNameOfHandle (CCTK_INT IN sym_handle)
REQUIRES FUNCTION SymmetryNameOfHandle
CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name)
REQUIRES FUNCTION SymmetryHandleOfName
CCTK_REAL FUNCTION GetScalProdCoeff ()
USES FUNCTION GetScalProdCoeff
SUBROUTINE Diff_gv ( CCTK_POINTER_TO_CONST IN cctkGH, \
CCTK_INT IN dir, \
CCTK_REAL IN ARRAY var, \
CCTK_REAL OUT ARRAY dvar, \
CCTK_INT IN table_handle )
USES FUNCTION Diff_gv
REAL pressuredensity TYPE=gf TAGS='tensortypealias="scalar"'
{
gappa
} "gamma=(rho+p)/p*dp/drho"
REAL tovsurface TYPE=scalar
{
TOV_surface
} "TOV_R+Xi(TOV_R)"
REAL metricpots TYPE=gf TAGS='tensortypealias="scalar"'
{
PHI, LAMBDA, drPHI, drLAMBDA, A, B, C, drrho, drpress
} "metric potentials: ds^2= -exp(2*PHI) dt^2 + exp(2*LAMBDA)dr^2 + ... and
depending quantities,"
REAL trafo TYPE=gf TAGS='tensortypealias="scalar"'
{
dxdr, dydr, dzdr, grid_r
} "coordiante transformation for derivatives... dxdrr, dydrr, dzdrr and
griddpoint radius"
REAL paramcheck TYPE=gf TAGS='tensortypealias="scalar"'
{
param_dr, deriv_error
} "to check accuracy derivatives"
REAL scalar TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
Xi
} "The scalar of the scalar wave equation fame"
REAL fluidscalar TYPE=gf
{
Xeta
} "fluid element displ."
REAL frequency TYPE=gf
{
omega, omegasquare
} "initial frequency"
REAL density TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
Pi
} "Time derivative of Xi"
REAL drscalar TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
drXi
} "Spatial derivatives of Xi"
REAL drdensity TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
drPi
} "Time derivative of Xi"
REAL scalardot TYPE=gf TAGS='tensortypealias="scalar" Prolongation="None"'
{
Xidot
} "RHS of Xi"
REAL densitydot TYPE=gf TAGS='tensortypealias="scalar" Prolongation="None"'
{
Pidot
} "RHS of Pi"
REAL velocity TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
H
} "H=drXi"
REAL drvelocity TYPE=gf TIMELEVELS=2 TAGS='tensortypealias="scalar"'
{
drH
} "Spatial derivatives of Xi"
REAL velocitydot TYPE=gf TAGS='tensortypealias="scalar" Prolongation="None"'
{
Hdot
} "RHS of H"
# Schedule definitions for thorn SWTNS
# Evolved variables
STORAGE: pressuredensity
STORAGE: metricpots
STORAGE: trafo
STORAGE: paramcheck
STORAGE: scalar[2] density[2] velocity[2]
STORAGE: drscalar[2] drvelocity[2]
STORAGE: drdensity[2]
STORAGE: tovradius
STORAGE: scalardot densitydot velocitydot
STORAGE: tovsurface
STORAGE: fluidscalar
STORAGE: frequency
#################################################
# Startup and initial
#################################################
SCHEDULE SWTNS_startup AT startup
{
LANG: C
OPTIONS: meta
} "Register banner with Cactus"
SCHEDULE SWTNS_register_MoL IN MoL_Register
{
LANG: C
OPTIONS: meta
} "Register variables with MoL"
SCHEDULE SWTNS_init AT initial AFTER TOV_C_Exact
{
LANG: Fortran
} "Initialise the system"
##################################################
# Calculate the RHS
##################################################
SCHEDULE SWTNS_calc_rhs IN MoL_CalcRHS
{
LANG: Fortran
} "Calculate the RHS"
#SCHEDULE SWTNS_surface IN MoL_CalcRHS AFTER SWTNS_calc_rhs
#{
# LANG: Fortran
#} "Calculate TOV radius"
##################################################
# Apply the boundary conditions
##################################################
SCHEDULE SWTNS_outerboundary IN MoL_PostStep
{
LANG: Fortran
} "Apply outer boundaries"
SCHEDULE SWTNS_RHS_outerboundary IN MoL_RHSBoundaries
{
LANG: Fortran
} "Apply MoL RHS outer boundaries (TOV_surface)"
SCHEDULE SWTNS_boundaries IN MoL_PostStep AFTER SWTNS_outerboundary
{
LANG: Fortran
OPTIONS: level
SYNC: scalar density drscalar scalardot densitydot velocity velocitydot
drvelocity drdensity
} "Select the boundary condition"
SCHEDULE GROUP ApplyBCs AS SWTNS_ApplyBCs IN MoL_PostStep AFTER SWTNS_boundaries
{
} "Apply boundary conditions"
##################################################
# Analysis: Frequencies, speed of sound?
##################################################
#if (recalculate_rhs)
#{
# SCHEDULE SWTNS_calc_rhs AT analysis
# {
# LANG: Fortran
# SYNC: scalardot densitydot
# STORAGE: scalardot densitydot
# TRIGGERS: scalardot densitydot
# } "Calculate the RHS"
#}
#todo
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
subroutine SWTNS_outerboundary (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: i, j, k, spec_i, spec_j, spec_k
do k = 1, cctk_lsh(3)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <= (TOV_surface + rprec) )then
H(i,j,k) = 0
!drXi(i,j,k) = 0
drPi(i,j,k) = 0
else if ( grid_r(i,j,k) > (TOV_surface + rprec) ) then
Xi(i,j,k) = 0
Pi(i,j,k) = 0
H(i,j,k) = 0
end if
end do
end do
end do
end subroutine SWTNS_outerboundary
subroutine SWTNS_RHS_outerboundary (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: i, j, k, spec_i, spec_j, spec_k
do k = 1, cctk_lsh(3)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
if (grid_r(i,j,k) > (TOV_surface + rprec) )then
Xidot(i,j,k) = 0
Hdot(i,j,k) = 0
Pidot(i,j,k) = 0
else if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <= (TOV_surface + rprec) )then
Hdot(i,j,k) = 0
! drXi(i,j,k) = 0
Pidot(i,j,k) = B(i,j,k)*drH(i,j,k) + C(i,j,k)*Xi(i,j,k)
end if
end do
end do
end do
end subroutine SWTNS_RHS_outerboundary
subroutine SWTNS_boundaries (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
character :: fbound*1000
integer :: fboundlen
integer :: ierr
!write(*,*) "SWTNS_boundaries"
call CCTK_FortranString (fboundlen, bound, fbound)
if (fboundlen > len(fbound)) call CCTK_WARN (0, "internal error")
ierr = Boundary_SelectGroupForBC &
(cctkGH, CCTK_ALL_FACES, +1, -1, "SWTNS::scalar", fbound)
if (ierr/=0) call CCTK_WARN (0, "internal error")
ierr = Boundary_SelectGroupForBC &
(cctkGH, CCTK_ALL_FACES, +1, -1, "SWTNS::density", fbound)
if (ierr/=0) call CCTK_WARN (0, "internal error")
ierr = Boundary_SelectGroupForBC &
(cctkGH, CCTK_ALL_FACES, +1, -1, "SWTNS::velocity", fbound)
if (ierr/=0) call CCTK_WARN (0, "internal error")
end subroutine SWTNS_boundaries
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"
subroutine SWTNS_init (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
integer :: i, j, k
!CCTK_REAL :: TOV_R global
CCTK_REAL :: pio, helper
pio = acos(-1.0d0)
TOV_surface = TOV_R
write(*,*)"Check if TOV_R is global..."
write(*,*)"Radius = ", TOV_R
write(*,*) "Surface = " , TOV_surface
write(*,*) "precission limit for grid points = " , rprec
! write(*,*) "grid boundary = " , gridboundary
write(*,*)"start initial"
!$OMP PARALLEL DO private(i,j,k)
do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1) !coordinate transformation factors
grid_r(i,j,k) = r(i,j,k)
if(patchnumber(i,j,k) > 0) then
dxdr(i,j,k) = iJ13(i,j,k)
dydr(i,j,k) = iJ23(i,j,k)
dzdr(i,j,k) = iJ33(i,j,k)
else if (patchnumber(i,j,k) == 0 .AND. grid_r(i,j,k) > 0) then
dxdr(i,j,k) = x(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
dydr(i,j,k) = y(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
dzdr(i,j,k) = z(i,j,k)/sqrt(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
end if
if(grid_r(i,j,k) == 0) then
dxdr(i,j,k) = 0
dydr(i,j,k) = 0
dzdr(i,j,k) = 0 !filling initial data starts here
PHI(i,j,k) = TOV_PHI(i,j,k)
drPHI(i,j,k) = 0
LAMBDA(i,j,k) = 0
drLAMBDA(i,j,k) = 0
drpress(i,j,k) = 0
drrho(i,j,k) = 0
A(i,j,k) = 0
B(i,j,k) = TOV_Gamma*press(i,j,k)/rho(i,j,k)*exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))
C(i,j,k) = 0
Xi(i,j,k) = (sin(pio*r(i,j,k)/TOV_surface)*cos(pio*r(i,j,k)/TOV_surface)-pio*r(i,j,k)/TOV_surface)*SWTNS_amplitude
Xeta(i,j,k) = 0
drXi(i,j,k) = pio * SWTNS_amplitude *(cos(2*pio*r(i,j,k)/TOV_surface)-1)/TOV_surface
H(i,j,k) = drXi(i,j,k)
drH(i,j,k) = -(2*pio**2*SWTNS_amplitude*sin(2*pio*r(i,j,k)/TOV_surface))/(TOV_surface**2)
omega(i,j,k) = 0
Pi(i,j,k) = 0
drPi(i,j,k) = 0
Xidot(i,j,k) = 0
Pidot(i,j,k) = 0
Hdot(i,j,k) = 0
else if (press(i,j,k) == 0 .OR. rho(i,j,k) == 0) then !outside the star
PHI(i,j,k) = TOV_PHI(i,j,k)
drPHI(i,j,k) = (TOV_mr(i,j,k) + 4*pio*r(i,j,k)**3 * press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
LAMBDA(i,j,k) = -0.5*log(1-2*TOV_mr(i,j,k)/r(i,j,k))
drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * rho(i,j,k)-TOV_mr(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
drpress(i,j,k) = -(rho(i,j,k) + press(i,j,k))*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3+press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
drrho(i,j,k) = 0 !manual approximation
A(i,j,k) = 0 !manual
B(i,j,k) = 0 !manual
C(i,j,k) = exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))*(drPHI(i,j,k)**2+4*drPHI(i,j,k)/r(i,j,k)-8*pio*exp(2*LAMBDA(i,j,k))*press(i,j,k) )
Xi(i,j,k) = 0
Xeta(i,j,k) = 0
drXi(i,j,k) = 0
H(i,j,k) = 0
drH(i,j,k) = 0
omega(i,j,k) = 0
Pi(i,j,k) = 0
drPi(i,j,k) = 0
Xidot(i,j,k) = 0
Pidot(i,j,k) = 0
Hdot(i,j,k) = 0
else
PHI(i,j,k) = TOV_PHI(i,j,k)
drPHI(i,j,k) = (TOV_mr(i,j,k) + 4*pio*r(i,j,k)**3 * press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
LAMBDA(i,j,k) = -0.5*log(1-2*TOV_mr(i,j,k)/r(i,j,k))
drLAMBDA(i,j,k) = (4*pio*r(i,j,k)**3 * rho(i,j,k)-TOV_mr(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
drpress(i,j,k) = -(rho(i,j,k) + press(i,j,k))*(TOV_mr(i,j,k)+4*pio*r(i,j,k)**3+press(i,j,k))/(r(i,j,k)*(r(i,j,k)-2*TOV_mr(i,j,k)))
drrho(i,j,k) = drpress(i,j,k)*rho(i,j,k)/(TOV_Gamma*press(i,j,k))
A(i,j,k) = (exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))/rho(i,j,k))*(TOV_Gamma*press(i,j,k)*(drLAMBDA(i,j,k) + 3*drPHI(i,j,k)) &
-(TOV_Gamma-1)*press(i,j,k)*drPHI(i,j,k)+TOV_Gamma*drpress(i,j,k) -2*TOV_Gamma*press(i,j,k)/r(i,j,k) )
B(i,j,k) = exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))*TOV_Gamma*press(i,j,k)/rho(i,j,k)
C(i,j,k) = exp(2*PHI(i,j,k)-2*LAMBDA(i,j,k))*(drPHI(i,j,k)**2+4*drPHI(i,j,k)/r(i,j,k)-8*pio*exp(2*LAMBDA(i,j,k))*press(i,j,k) )
Xi(i,j,k) = (sin(pio*r(i,j,k)/TOV_surface)*cos(pio*r(i,j,k)/TOV_surface)-pio*r(i,j,k)/TOV_surface)*SWTNS_amplitude
Xeta(i,j,k) = Xi(i,j,k)*exp(LAMBDA(i,j,k)+PHI(i,j,k))/(r(i,j,k)**2)
drXi(i,j,k) = pio * SWTNS_amplitude *(cos(2*pio*r(i,j,k)/TOV_surface)-1)/TOV_surface
H(i,j,k) = drXi(i,j,k)
drH(i,j,k) = -(2*pio**2*SWTNS_amplitude*sin(2*pio*r(i,j,k)/TOV_surface))/(TOV_surface**2)
omegasquare(i,j,k) = 0!(-1)*(A(i,j,k)*drXi(i,j,k)/Xi(i,j,k) + B(i,j,k)*drH(i,j,k)/Xi(i,j,k)+C(i,j,k))
if(omegasquare(i,j,k) <= 0) then
omega(i,j,k) = 0
Pi(i,j,k) = 0
drPi(i,j,k) = 0
Xidot(i,j,k) = 0
Pidot(i,j,k) = 0
Hdot(i,j,k) = 0
else
omega(i,j,k) = sqrt(omegasquare(i,j,k))
Pi(i,j,k) = Xi(i,j,k)*omega(i,j,k)
drPi(i,j,k) = drXi(i,j,k)*omega(i,j,k)
Xidot(i,j,k) = Xi(i,j,k)*omega(i,j,K)
Hdot(i,j,k) = drXi(i,j,k)*omega(i,j,k)
Pidot(i,j,k) = -Xi(i,j,k)*omegasquare(i,j,k)
end if
end if
end do
end do
end do
!$OMP END PARALLEL DO
write(*,*) "Coordinate transformation factors for spatial derivatives have been prepared for Thornburg04 and Thornburg04nc."
write(*,*) "Calculated metric potentials, derivatives and coefficients A,B,C"
write(*,*)"done initial!"
call SWTNS_drparainit(CCTK_ARGUMENTS)
end subroutine SWTNS_init
subroutine SWTNS_drparainit (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
CCTK_INT, parameter :: izero = 0
! CCTK_INT :: i, j, k
integer, parameter :: ik = kind (izero)
integer :: na, nb, nc
CCTK_REAL, dimension(:,:,:), allocatable &
:: param_dx, param_dy, param_dz
na = cctk_lsh(1); nb = cctk_lsh(2); nc = cctk_lsh(3)
allocate( param_dx(na,nb,nc), &
param_dy(na,nb,nc), &
param_dz(na,nb,nc) )
!check derivative algorithem with drPHI
call globalDiff_gv (cctkGH, 0_ik, PHI, param_dx, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
call globalDiff_gv (cctkGH, 1_ik, PHI, param_dy, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
call globalDiff_gv (cctkGH, 2_ik, PHI, param_dz, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
param_dr = dxdr * param_dx + dydr * param_dy + dzdr * param_dz
deriv_error = drPHI - param_dr
deallocate (param_dx, param_dy,param_dz)
!write(*,*) "Calculated drP"
end subroutine SWTNS_drparainit
#include <assert.h>
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
int
SWTNS_startup (void)
{
int ierr;
ierr = CCTK_RegisterBanner ("SWTNS");
assert (! ierr);
return 0;
}
void
SWTNS_register_MoL (CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS;
DECLARE_CCTK_PARAMETERS;
int evolved_group, rhs_group;
int ierr;
evolved_group = CCTK_GroupIndex ("SWTNS::scalar");
assert (evolved_group >= 0);
rhs_group = CCTK_GroupIndex ("SWTNS::scalardot");
assert (rhs_group >= 0);
ierr = MoLRegisterEvolvedGroup (evolved_group, rhs_group);
assert (! ierr);
evolved_group = CCTK_GroupIndex ("SWTNS::velocity");
assert (evolved_group >= 0);
rhs_group = CCTK_GroupIndex ("SWTNS::velocitydot");
assert (rhs_group >= 0);
ierr = MoLRegisterEvolvedGroup (evolved_group, rhs_group);
assert (! ierr);
evolved_group = CCTK_GroupIndex ("SWTNS::density");
assert (evolved_group >= 0);
rhs_group = CCTK_GroupIndex ("SWTNS::densitydot");
assert (rhs_group >= 0);
ierr = MoLRegisterEvolvedGroup (evolved_group, rhs_group);
assert (! ierr);
}
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"
subroutine SWTNS_calc_rhs (CCTK_ARGUMENTS)
implicit none
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_FUNCTIONS
DECLARE_CCTK_PARAMETERS
CCTK_INT, parameter :: izero = 0
CCTK_INT :: i, j, k, spec_i, spec_j, spec_k
CCTK_REAL :: pio
integer, parameter :: ik = kind (izero)
integer :: na, nb, nc
CCTK_REAL, dimension(:,:,:), allocatable &
:: TOVold, dxXi, dyXi, dzXi, &
dxPi, dyPi, dzPi, dxH, dyH, dzH
na = cctk_lsh(1); nb = cctk_lsh(2); nc = cctk_lsh(3)
! write(*,*)"allocate RHS..."
pio=acos(-1.0d0)
!write(*,*) "Surface = " , TOV_surface
allocate( dxXi(na, nb, nc), &
dyXi(na, nb, nc), &
dzXi(na, nb, nc), &
dxPi(na, nb, nc), &
dyPi(na, nb, nc), &
dzPi(na, nb, nc), &
dxH(na, nb, nc), &
dyH(na, nb, nc), &
dzH(na, nb, nc) )
!call SWTNS_derivatives (CCTK_ARGUMENTS)
!!! Wave Eq.
! Xidot = Pi
! Hdot = drPi
! Pidot = A*H + B*drH + C*Xi
!write(*,*) "Define Xeta"
!Boundary Xi(r=0) = 0
do k = 1, cctk_lsh(3)
do j = 1, cctk_lsh(2)
do i = 1, cctk_lsh(1)
if (grid_r(i,j,k) == 0 ) then
Xi(i,j,k) = 0
Xeta(i,j,k) = 0
else
Xeta(i,j,k) = Xi(i,j,k)*exp(PHI(i,j,k))/(r(i,j,k)**2)
end if
end do
end do
end do
!updating TOV surface in PostStep
!Updating TOV_surface
surfaceloop: do k=1,cctk_lsh(3)
do j=1,cctk_lsh(2)
do i=1,cctk_lsh(1)
if(patchnumber(i,j,k) == 1) then
if (grid_r(i,j,k) >= (TOV_surface - rprec) .AND. grid_r(i,j,k) <= (TOV_surface) .AND. Xi(i,j,k) /= 0) then
TOV_surface = TOV_R + Xeta(i,j,k) !or Xi
exit surfaceloop
end if
else
exit surfaceloop
end if
end do
end do
end do surfaceloop
!Ableitung drXi
! call globalDiff_gv (cctkGH, 0_ik, Xi, dxXi, J11, J21, J31, &
! J12, J22, J32, J13, J23, J33,-1_ik)
! call globalDiff_gv (cctkGH, 1_ik, Xi, dyXi, J11, J21, J31, &
! J12, J22, J32, J13, J23, J33, -1_ik)
! call globalDiff_gv (cctkGH, 2_ik, Xi, dzXi, J11, J21, J31, &
! J12, J22, J32, J13, J23, J33, -1_ik)
! drXi = dxXi*dxdr + dyXi*dydr + dzXi*dzdr
!Ableitung drPi
call globalDiff_gv (cctkGH, 0_ik, Pi, dxPi, J11, J21, J31, &
J12, J22, J32, J13, J23, J33,-1_ik)
call globalDiff_gv (cctkGH, 1_ik, Pi, dyPi, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
call globalDiff_gv (cctkGH, 2_ik, Pi, dzPi, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
drPi = dxPi*dxdr + dyPi*dydr + dzPi*dzdr
!Ableitung drH
call globalDiff_gv (cctkGH, 0_ik, H, dxH, J11, J21, J31, &
J12, J22, J32, J13, J23, J33,-1_ik)
call globalDiff_gv (cctkGH, 1_ik, H, dyH, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
call globalDiff_gv (cctkGH, 2_ik, H, dzH, J11, J21, J31, &
J12, J22, J32, J13, J23, J33, -1_ik)
drH = dxH*dxdr + dyH*dydr + dzH*dzdr
!!! Wave Eq.
Pidot = A*H + B*drH + C*Xi
Hdot = drPi
Xidot = Pi
deallocate ( dxXi, dyXi, dzXi, dxPi, dyPi, dzPi, dxH, dyH, dzH)
end subroutine SWTNS_calc_rhs
_______________________________________________
Users mailing list
[email protected]
http://lists.einsteintoolkit.org/mailman/listinfo/users
