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