# [Users] Boundary issues solving wave equations with MoL

Severin Frank severin.frank at uni-tuebingen.de
Sat Dec 21 14:39:42 CST 2019

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

Ẍ=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

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.einsteintoolkit.org/pipermail/users/attachments/20191221/04a98185/attachment-0001.html
-------------- next part --------------
#===============================================================
# Testing SWTNS on a static TOV star
#
# Severin Frank
#===============================================================

#################### Thorns ####################################
ActiveThorns = "Time MoL"

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

SWTNS::rprec 				= 0.1
TOVSolver::tovrprec 	                = 0.1
SWTNS::SWTNS_amplitude			= 0.001 #0.0031830989

Coordinates::n_angular                  = 9

Driver::ghost_size                      = 3
Coordinates::patch_boundary_size        = 3

Coordinates::outer_boundary_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

TOVSolver::TOV_Rho_Central = 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

#################### 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}
"
-------------- 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_dbdx, \
CCTK_REAL IN ARRAY J_dcdx, \
CCTK_REAL IN ARRAY J_dbdy, \
CCTK_REAL IN ARRAY J_dcdy, \
CCTK_REAL IN ARRAY J_dbdz, \
CCTK_REAL IN ARRAY J_dcdz, \
CCTK_REAL IN ARRAY dJ_dbdxdx, \
CCTK_REAL IN ARRAY dJ_dcdxdx, \
CCTK_REAL IN ARRAY dJ_dbdxdy, \
CCTK_REAL IN ARRAY dJ_dcdxdy, \
CCTK_REAL IN ARRAY dJ_dbdxdz, \
CCTK_REAL IN ARRAY dJ_dcdxdz, \
CCTK_REAL IN ARRAY dJ_dbdydy, \
CCTK_REAL IN ARRAY dJ_dcdydy, \
CCTK_REAL IN ARRAY dJ_dbdydz, \
CCTK_REAL IN ARRAY dJ_dcdydz, \
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"
-------------- next part --------------
# Schedule definitions for thorn SWTNS

# Evolved variables
STORAGE: pressuredensity
STORAGE: metricpots
STORAGE: trafo
STORAGE: paramcheck
STORAGE: scalar density velocity
STORAGE: drscalar drvelocity
STORAGE: drdensity
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

##################################################
# 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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: boundary.F90
Type: text/x-fortran
Size: 2307 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191221/04a98185/attachment-0004.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: init.F90
Type: text/x-fortran
Size: 6528 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191221/04a98185/attachment-0005.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: register.c
Type: text/x-csrc
Size: 1107 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191221/04a98185/attachment-0006.bin
-------------- next part --------------
A non-text attachment was scrubbed...
Name: rhs.F90
Type: text/x-fortran
Size: 3518 bytes
Desc: not available
Url : http://lists.einsteintoolkit.org/pipermail/users/attachments/20191221/04a98185/attachment-0007.bin