<div dir="ltr">Hi Roland, thank you so much for bringing this to the ET call.<div><br></div><div>I see. OK, I&#39;ll ask Oleg.</div><div><br></div><div>Thanks!<br clear="all"><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><font color="#999999">Luciano</font></div></div></div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Mar 9, 2023 at 12:37 PM Roland Haas &lt;<a href="mailto:rhaas@illinois.edu">rhaas@illinois.edu</a>&gt; wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Hello Luciano,<br>
<br>
We quickly discussed this in last week&#39;s ET call.<br>
<br>
Unfortunately not one present had any experience using Llaman with<br>
GRHydro and a non-trivial topology.<br>
<br>
Oleg Korobkin did some similar work on accretion disks using a Cactus<br>
based code with a central cutout (but not Llama based but using an<br>
earlier version of Cactus/Carpet based multiblock): <br>
<br>
<a href="https://arxiv.org/abs/1011.3010" rel="noreferrer" target="_blank">https://arxiv.org/abs/1011.3010</a><br>
<br>
You may be a able to reach out to him (at LANL).<br>
<br>
My own recollection is that hydro_excision is not used very frequently<br>
in GRHydro, to it may or may not actually work at all.<br>
<br>
Note when you say &quot;origin&quot; since this is the Llama grid with a hole in<br>
the center, you mean &quot;inner cutout&quot; (ie not r=0 but instead r=3.0)?<br>
<br>
Yours,<br>
Roland<br>
<br>
&gt; Hi everyone,<br>
&gt; <br>
&gt; I started testing the Thornburg04nc coordinates (spherical patches without<br>
&gt; Cartesian box) from Llama with (a modified version) of GRHydro. I started<br>
&gt; with just an atmosphere around a Kerr-Schild metric evolving only hydro.<br>
&gt; I&#39;m seeing instabilities near the origin, so I was wondering if someone<br>
&gt; already encountered this issue. In particular, I wanted to know if someone<br>
&gt; has a good idea of how to treat the inner boundary. Right now, I&#39;m putting<br>
&gt; the inner boundary inside the horizon, above the ring singularity, and<br>
&gt; applying hydro excision. I&#39;ve tried a combination of things playing with<br>
&gt; those parameters but without success. Has anyone tried a similar thing<br>
&gt; evolving hydro? Just want to double check if there&#39;s something obvious in<br>
&gt; the grid setup I&#39;m missing.<br>
&gt; <br>
&gt; Below I attach the relevant part of the par file I&#39;m using. Main difference<br>
&gt; with the public version of GRHydro is that I&#39;m using Siegel+2017 c2p which<br>
&gt; works just fine with the same setup in Cartesian coordinates.<br>
&gt; <br>
&gt; Thanks a lot!!<br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; ######<br>
&gt; <br>
&gt; ActiveThorns = &quot;Time MoL&quot;<br>
&gt; ActiveThorns = &quot;Coordbase CartGrid3d Boundary StaticConformal&quot;<br>
&gt; ActiveThorns = &quot;SymBase ADMBase TmunuBase HydroBase InitBase ADMCoupling<br>
&gt; ADMMacros&quot;<br>
&gt; ActiveThorns = &quot;IOUtil Formaline&quot;<br>
&gt; ActiveThorns = &quot;SpaceMask CoordGauge Constants LocalReduce aeilocalinterp<br>
&gt; LoopControl&quot;<br>
&gt; ActiveThorns = &quot;Carpet CarpetLib CarpetReduce CarpetRegrid2 &quot;<br>
&gt; ActiveThorns = &quot;Coordinates GenericFD Interpolate2 CarpetInterp2<br>
&gt; CarpetInterp&quot;<br>
&gt; ActiveThorns = &quot;CarpetIOASCII CarpetIOScalar CarpetIOHDF5 CarpetIOBasic&quot;<br>
&gt; <br>
&gt; ################################################################################<br>
&gt; # Grid setup and general parameters<br>
&gt; ################################################################################<br>
&gt; <br>
&gt; # Grid<br>
&gt; <br>
&gt; CartGrid3D::type               = &quot;multipatch&quot;<br>
&gt; Coordinates::coordinate_system = &quot;Thornburg04nc&quot;<br>
&gt; CartGrid3D::set_coordinate_ranges_on    = &quot;all maps&quot;<br>
&gt; <br>
&gt; Coordinates::h_radial            =   0.35<br>
&gt; Coordinates::sphere_inner_radius =  3.0<br>
&gt; Coordinates::sphere_outer_radius =  90.0<br>
&gt; Coordinates::n_angular           =  25<br>
&gt; <br>
&gt; #Coordinates::radial_stretch             = yes<br>
&gt; #Coordinates::stretch_rmin_1             = 10.0<br>
&gt; #Coordinates::stretch_rmax_1             = 90.0<br>
&gt; <br>
&gt; Coordinates::stagger_patch_boundaries = &quot;no&quot;<br>
&gt; Coordinates::stagger_outer_boundaries = &quot;no&quot;<br>
&gt; <br>
&gt; Interpolate::interpolator_order = 4<br>
&gt; Interpolate2::shift_edges        = no<br>
&gt; <br>
&gt; Coordinates::store_jacobian            = yes # I need to<br>
&gt; Coordinates::store_inverse_jacobian    = yes # Check for runtime errors<br>
&gt; #Coordinates::store_jacobian_derivative = yes<br>
&gt; Coordinates::store_volume_form         = yes<br>
&gt; <br>
&gt; #GenericFD::jacobian_group            = &quot;Coordinates::jacobian&quot;<br>
&gt; #GenericFD::jacobian_derivative_group = &quot;Coordinates::jacobian2&quot;<br>
&gt; <br>
&gt; # General Carpet parameters:<br>
&gt; Carpet::enable_all_storage       = &quot;no&quot;<br>
&gt; Carpet::use_buffer_zones         = &quot;yes&quot;<br>
&gt; Carpet::use_overlap_zones        = &quot;no&quot;<br>
&gt; Carpet::schedule_barriers        = &quot;no&quot;<br>
&gt; #Carpet::processor_topology       = &quot;recursive&quot; # recursive breaks<br>
&gt; multipatch<br>
&gt; <br>
&gt; Carpet::check_for_poison         = &quot;no&quot;<br>
&gt; Carpet::init_3_timelevels        = &quot;no&quot;<br>
&gt; Carpet::init_fill_timelevels     = &quot;yes&quot;<br>
&gt; <br>
&gt; CarpetLib::poison_new_memory         = &quot;yes&quot;<br>
&gt; CarpetLib::poison_value              = 114<br>
&gt; CarpetLib::check_bboxes              = &quot;no&quot;<br>
&gt; CarpetLib::interleave_communications = &quot;yes&quot;<br>
&gt; CarpetLib::combine_sends             = &quot;no&quot; # Erik says this is faster and<br>
&gt; does not use more memory than yes<br>
&gt; CarpetLib::combine_recompose         = &quot;no&quot; # yes is the default, set to no<br>
&gt; if we run out of memory during recomposing<br>
&gt; #CarpetLib::restriction_order_space   = 0<br>
&gt; <br>
&gt; <br>
&gt; CarpetInterp::tree_search = &quot;yes&quot;<br>
&gt; CarpetInterp::check_tree_search = &quot;no&quot;<br>
&gt; <br>
&gt; # OVERLAP<br>
&gt; Coordinates::additional_overlap_size = 3<br>
&gt; <br>
&gt; # BOUNDARY ZONES (needs to be set to three)<br>
&gt; Coordinates::patch_boundary_size = 3<br>
&gt; Coordinates::outer_boundary_size = 3<br>
&gt; #Carpet::prolongation_order_space = 5<br>
&gt; #Carpet::prolongation_order_time  = 2<br>
&gt; Carpet::refinement_centering     = &quot;vertex&quot;<br>
&gt; <br>
&gt; #====================================<br>
&gt; # TmunuBase &amp; HydroBase<br>
&gt; #====================================<br>
&gt; TmunuBase::stress_energy_storage = yes<br>
&gt; TmunuBase::stress_energy_at_RHS  = yes<br>
&gt; TmunuBase::timelevels            =  1<br>
&gt; TmunuBase::prolongation_type     = none<br>
&gt; HydroBase::timelevels            = 3<br>
&gt; <br>
&gt; #############################################################<br>
&gt; # Carpet<br>
&gt; #############################################################<br>
&gt; <br>
&gt; Carpet::domain_from_multipatch = yes<br>
&gt; Carpet::ghost_size             = 3<br>
&gt; # not for FD<br>
&gt; #Carpet::granularity            = 5   # (order+1)*m<br>
&gt; #Carpet::granularity_boundary   = 1<br>
&gt; #Carpet::poison_new_timelevels = &quot;yes&quot;   # does not work with persistent RHS<br>
&gt; <br>
&gt; #--------CarpetMask--------------------------<br>
&gt; ActiveThorns = &quot;CarpetMask&quot;<br>
&gt; CarpetMask::excluded_surface       [0] = 1<br>
&gt; CarpetMask::excluded_surface_factor[0] = 1.0<br>
&gt; #<br>
&gt; #############################################################<br>
&gt; # Time integration<br>
&gt; #############################################################<br>
&gt; <br>
&gt; Cactus::terminate       = &quot;time&quot;<br>
&gt; Cactus::cctk_final_time = 2000000000000<br>
&gt; #Cactus::cctk_itlast = 1<br>
&gt; <br>
&gt; Time::timestep_method = &quot;given&quot;<br>
&gt; Time::timestep        = 0.10<br>
&gt; <br>
&gt; MethodOfLines::ODE_method             = &quot;RK4&quot;<br>
&gt; MethodOfLines::MoL_Intermediate_Steps = 4<br>
&gt; MethodOfLines::MoL_Num_Scratch_Levels = 1<br>
&gt; MethodOfLines::MoL_NaN_Check = no<br>
&gt; <br>
&gt; <br>
&gt; ################################################################################<br>
&gt; # Hydro evolution<br>
&gt; ################################################################################<br>
&gt; <br>
&gt; SpaceMask::use_mask      = &quot;yes&quot;<br>
&gt; <br>
&gt; #====================================<br>
&gt; # GRHydro<br>
&gt; #====================================<br>
&gt; ActiveThorns = &quot;EOS_Omni&quot;<br>
&gt; ActiveThorns = &quot;GRHydro&quot;<br>
&gt; <br>
&gt; HydroBase::evolution_method      = &quot;GRHydro&quot;<br>
&gt; HydroBase::Bvec_evolution_method = &quot;GRHydro&quot;<br>
&gt; # use c++ version of GRHydro<br>
&gt; GRHydro::use_cxx_code = &quot;yes&quot;<br>
&gt; <br>
&gt; GRHydro::method_type               = &quot;RSA FV&quot;<br>
&gt; GRHydro::riemann_solver            = &quot;HLLE&quot;<br>
&gt; GRHydro::recon_method              = &quot;ppm&quot;<br>
&gt; GRHydro::recon_vars                = &quot;primitive&quot;<br>
&gt; GRHydro::GRHydro_stencil           = 3<br>
&gt; GRHydro::bound                     = &quot;flat&quot; #also tried none<br>
&gt; GRHydro::rho_abs_min               = 1e-15<br>
&gt; GRHydro::initial_rho_abs_min       = 5e-14<br>
&gt; GRHydro::GRHydro_atmo_tolerance    = 0.01<br>
&gt; <br>
&gt; <br>
&gt; # -- settings for 3-param (tabulated) EOS --<br>
&gt; grhydro::grhydro_eos_type                =  &quot;General&quot;<br>
&gt; grhydro::grhydro_eos_table               =  &quot;nuc_eos&quot;<br>
&gt; EOS_Omni::nuceos_read_table = &quot;yes&quot;<br>
&gt; EOS_Omni::nuceos_table_name =<br>
&gt; &quot;/home/d/dsiegel/lcombi/repository/EOS_tables/etk_bns_evolution/sfho/Hempel_SFHoEOS_rho222_temp180_ye60_version_1.1_20120817.h5&quot;<br>
&gt; GRHydro::GRHydro_eos_rf_prec = 1e-10<br>
&gt; hydrobase::Y_e_evolution_method          =  &quot;GRHydro&quot;<br>
&gt; hydrobase::temperature_evolution_method  =  &quot;GRHydro&quot;<br>
&gt; <br>
&gt; GRHydro::reconstruct_temper             =  &quot;yes&quot; # more robust for<br>
&gt; reconstruction<br>
&gt; GRHydro::GRHydro_hot_atmo_temp          =  0.012 #8.62e-06=1e5K<br>
&gt; GRHydro::GRHydro_hot_atmo_Y_e           =  0.5<br>
&gt; grhydro::grhydro_y_e_min               =  0.012<br>
&gt; grhydro::grhydro_y_e_max               =  0.599<br>
&gt; <br>
&gt; <br>
&gt; GRHydro::Grhydro_MaxNumConstrainedVars = 34 #33<br>
&gt; GRHydro::GRHydro_MaxNumEvolvedVars     = 10<br>
&gt; <br>
&gt; # -- settings for con2prim --<br>
&gt; GRHydro::GRHydro_c2pMHD_method = 1<br>
&gt; GRHydro::GRHydro_c2pMHD_method_backup = 5<br>
&gt; GRHydro::GRHydro_use_c2pMHD_method_backup = &quot;yes&quot;<br>
&gt; GRHydro::GRHydro_c2p_grace_radius  = 7.5<br>
&gt; grhydro::grhydro_c2p_warnlevel          = 1 #1<br>
&gt; grhydro::grhydro_c2p_warn_from_reflevel = 5 #5 # 6<br>
&gt; <br>
&gt; # --- EOS Omni ---<br>
&gt; eos_omni::poly_gamma                            = 1.33333333<br>
&gt; #eos_omni::poly_gamma_initial                    = 1.333333333<br>
&gt; eos_omni::poly_k                                = 0.3605 #145.877704 #100<br>
&gt; eos_omni::gl_gamma                              = 1.33333333<br>
&gt; eos_omni::gl_k                                  = 0.3605 #145.877704 #100<br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; ################################################################################<br>
&gt; # Spacetime evolution<br>
&gt; ################################################################################<br>
&gt; <br>
&gt; ADMBase::metric_type     = &quot;physical&quot;<br>
&gt; ADMBase::evolution_method        = &quot;static&quot;<br>
&gt; ADMBase::lapse_evolution_method  = &quot;static&quot;<br>
&gt; ADMBase::shift_evolution_method  = &quot;static&quot;<br>
&gt; ADMBase::dtlapse_evolution_method= &quot;static&quot;<br>
&gt; ADMBase::dtshift_evolution_method= &quot;static&quot;<br>
&gt; <br>
&gt; ADMBase::metric_timelevels = 3<br>
&gt; ADMBase::lapse_timelevels = 3<br>
&gt; ADMBase::shift_timelevels = 3<br>
&gt; ADMMacros::spatial_order = 4<br>
&gt; <br>
&gt; #====================================<br>
&gt; # Apparent horizon &amp; excision<br>
&gt; #====================================<br>
&gt; ActiveThorns = &quot;hydro_initexcision&quot;<br>
&gt; hydro_initexcision::hydro_initexcision  = &quot;yes&quot;<br>
&gt; hydro_initexcision::hydro_initexcision_type  = &quot;sphere&quot;<br>
&gt; hydro_initexcision::hydro_initexcision_coordinate_length  = 4.2 # should be<br>
&gt; slightly larger than a_BH = spin_BH*mass_BH<br>
&gt; <br>
&gt; # --- Settings for apparent horizon ---<br>
&gt; ActiveThorns = &quot;AHFinderDirect SetMask_SphericalSurface&quot;<br>
&gt; AHFinderDirect::N_horizons                               = 1<br>
&gt; AHFinderDirect::find_after_individual_time[1]            = 0<br>
&gt; AHFinderDirect::find_every                               = 20000<br>
&gt; AHFinderDirect::output_h_every                           = 20000<br>
&gt; AHFinderDirect::initial_guess_method[1]                   =<br>
&gt; &quot;Kerr/Kerr-Schild&quot;<br>
&gt; AHFinderDirect::initial_guess__Kerr_KerrSchild__x_posn[1] = 0.0<br>
&gt; AHFinderDirect::initial_guess__Kerr_KerrSchild__y_posn[1] = 0.0<br>
&gt; AHFinderDirect::initial_guess__Kerr_KerrSchild__z_posn[1] = 0.0<br>
&gt; AHFinderDirect::initial_guess__Kerr_KerrSchild__mass[1] = 3.0  # adjust to<br>
&gt; settings for TorusID<br>
&gt; AHFinderDirect::initial_guess__Kerr_KerrSchild__spin[1] = 0.8  # adjust to<br>
&gt; settings for TorusID<br>
&gt; AHFinderDirect::geometry_interpolator_name               = &quot;Lagrange<br>
&gt; polynomial interpolation&quot;<br>
&gt; AHFinderDirect::geometry_interpolator_pars               = &quot;order=4&quot;<br>
&gt; AHFinderDirect::surface_interpolator_name                = &quot;Lagrange<br>
&gt; polynomial interpolation&quot;<br>
&gt; AHFinderDirect::surface_interpolator_pars                = &quot;order=4&quot;<br>
&gt; AHFinderDirect::mask_buffer_thickness                    = 0<br>
&gt; AHFinderDirect::mask_radius_multiplier                   = 0.8<br>
&gt; AHFinderDirect::mask_radius_offset                       = 0<br>
&gt; AHFinderDirect::max_allowable_horizon_radius[1]          = 10<br>
&gt; AHFinderDirect::max_allowable_Theta_growth_iterations    = 10<br>
&gt; AHFinderDirect::max_allowable_Theta_nonshrink_iterations = 10<br>
&gt; AHFinderDirect::max_Newton_iterations__initial           = 50<br>
&gt; AHFinderDirect::max_Newton_iterations__subsequent        = 50<br>
&gt; AHFinderDirect::which_surface_to_store_info_by_name  [1] = &quot;apparent<br>
&gt; horizon&quot;<br>
&gt; <br>
&gt; ################################################################################<br>
&gt; # Initial data<br>
&gt; ################################################################################<br>
&gt; <br>
&gt; InitBase::initial_data_setup_method = &quot;init_some_levels&quot;<br>
&gt; ADMBase::initial_data            = &quot;exact&quot; #TorusID&quot;<br>
&gt; ADMBase::initial_lapse           = &quot;exact&quot; #TorusID&quot;<br>
&gt; ADMBase::initial_shift           = &quot;exact&quot; #TorusID&quot;<br>
&gt; ADMBase::initial_dtlapse         = &quot;zero&quot;<br>
&gt; ADMBase::initial_dtshift         = &quot;zero&quot;<br>
&gt; <br>
&gt; ActiveThorns = &quot;Exact&quot;<br>
&gt; Exact::exact_model = &quot;Kerr/Kerr-Schild&quot;<br>
&gt; Exact::Kerr_KerrSchild__mass     = 3.0<br>
&gt; Exact::Kerr_KerrSchild__parabolic = &quot;yes&quot;<br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; <br>
&gt; ---<br>
&gt; *Dr. Luciano Combi*<br>
&gt; Postdoctoral Researcher<br>
&gt; Perimeter Institute for Theoretical Physics<br>
&gt; CITA National Fellow (U. of Guelph)<br>
&gt; ---<br>
<br>
-- <br>
My email is as private as my paper mail. I therefore support encrypting<br>
and signing email messages. Get my PGP key from <a href="http://pgp.mit.edu" rel="noreferrer" target="_blank">http://pgp.mit.edu</a> .<br>
</blockquote></div>