[Users] GlobalDerivative becomes 0 for certain grid points (Llama)

Severin Frank severin.frank at uni-tuebingen.de
Wed Dec 4 11:13:42 CST 2019


Thank you Erik for your last reply and your helpful advice!

In the last days I changed a lot of my code and I noticed another 
(different?) issue:

I'm calculating spatial derivatives in r-direction for different grid 
functions by using Llama GlobalDerivative thorn, such as this example:

              call globalDiff_gv (cctkGH, 0, PHI, param_dx, J11, J21, J31, &
                                   J12, J22, J32, J13, J23, J33, -1_ik)
              call globalDiff_gv (cctkGH, 1, PHI, param_dy, J11, J21, J31, &
                                   J12, J22, J32, J13, J23, J33, -1_ik)
              call globalDiff_gv (cctkGH, 2, 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

, where dxdr, dydr, dzdr are the specific coordinate transformations and 
my coordinates are again Thornburg04 multipatch.

I noticed that I receive very odd (and wrong) zero-values on specific 
grid points in the positive x-sphere. Not only for param_dr, but also 
for param_dx, param_dy as well as param_dz. I should mention, that I did 
not change the GlobalDerivative thorn and this problem also appears, 
when calculating the derivatives of other variables then PHI. Also PHI 
itself is continuous in this area. I also initiate all "param"-variables 
with 1 so that the call of GlobalDerivative is actively changing them 
(to 0).

The same error also appears, when calling GlobalDiff directly with the 
inverse Jacobians like:

call globalDiff_gv (cctkGH, 2, PHI, param_dr, J11, J21, J31, J12, J22, 
J32, iJ13, iJ23, iJ33, -1_ik).

Furthermore it makes no difference whether this is scheduled at initial 
or in MoL_CalcRHS.

I attached plots from PHI, param_dx and param_dr so that you can see 
where those points are on the x-z plane. These are for an angular 
resolution of Coordinates::n_angular = 9. If I increase the angular 
resolution this affects also more grid points along the same radius (for 
example 8 points in the x-z plane for n_angular = 15).

Since I'm assuming, that the GobalDerivative thorn is working correct I 
don't know where this issue might come from. Does anyone have an idea 
for this?

Thank you all a lot!

Best regards,
Severin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: init.F90
Type: text/x-fortran
Size: 7122 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191204/83eb9558/attachment-0002.bin 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rhs.F90
Type: text/x-fortran
Size: 5516 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191204/83eb9558/attachment-0003.bin 
-------------- next part --------------
# 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
#################################################
# 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"

#SCHEDULE SWTNS_inittest IN CCTK_PREREGRID
#{
#  LANG: Fortran
#} "TEST the system"

 
#SCHEDULE SWTNS_regridtest AT CCTK_POSTREGRID AFTER SpatialCoordinates
#{
#  LANG: Fortran
#} "TEST 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 velocity drscalar drvelocity scalardot drdensity densitydot velocitydot
#} "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
-------------- next part --------------
# 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, param_dx, param_dy, param_dz
} "to check critical 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 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"
-------------- next part --------------
#===============================================================
# 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 CarpetReduce"

ActiveThorns = "IOUtil CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic"

ActiveThorns = "NaNChecker"

ActiveThorns = "Coordinates CoordinatesSymmetry Interpolate2 CarpetInterp2 GlobalDerivative SWTNS" 


#################### 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

SWTNS::gridboundary                     = 2 # ueberfluessig
Coordinates::sphere_inner_radius        = 1 #
Coordinates::sphere_outer_radius        = 10	#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


SWTNS::outer_bound ="dirichlet"

#################### 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::H
			  SWTNS::Hdot SWTNS::drH SWTNS::drPi
			  SWTNS::Xi SWTNS::Xidot SWTNS::Pi SWTNS::Pidot 
			  SWTNS::Xeta SWTNS::drXi SWTNS::param_dr SWTNS::param_dx SWTNS::param_dy SWTNS::param_dz"

#################### 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::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::Seve_in_cube{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::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::Pidot
					  SWTNS::drXi
					  SWTNS::drH
				          SWTNS::H
					  SWTNS::Hdot
				          SWTNS::grid_r{out_every=1000000000000000}
					  SWTNS::param_dr
					  SWTNS::param_dx{out_every=1000000000000000}
					  SWTNS::param_dy{out_every=1000000000000000}
					  SWTNS::param_dz{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}                             
                                           "                                          
-------------- next part --------------
A non-text attachment was scrubbed...
Name: param_dr.jpg
Type: image/jpeg
Size: 298079 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191204/83eb9558/attachment-0003.jpg 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: param_dx.jpg
Type: image/jpeg
Size: 320152 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191204/83eb9558/attachment-0004.jpg 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: PHI.jpg
Type: image/jpeg
Size: 298791 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191204/83eb9558/attachment-0005.jpg 


More information about the Users mailing list