<div dir="ltr"><div>Roland,</div><div> thank you for your precious suggestions.</div><div><br></div><div>I managed to solve the issue. After setting AHFinderDirect verbose level to "algorithm debug" I found out that some points of the horizon surface fell on an excised region:</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>AHFinderDirect::AHFinderDirect_find_horizons<br> in thorn AHFinderDirect, file /marconi/home/userexternal/fcattori/EinsteinToolkit/ET_2018_09/Cactus/arrangements/EinsteinAnalysis/AHFinderDirect/src/gr/expansion.cc:949:<br> -> interpolate_geometry():<br> one or more points on the trial horizon surface point<br> is/are in an excised region (or too close to the excision boundary)</div></blockquote><div><br></div><div>Then I added to the .par file the string</div><div><br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>AHFinderDirect::mask_is_noshrink = "false"<br></div></blockquote><div><br></div><div>and now everything is working properly: the surface is valid (sf_valid[2] = 1 from t_merger onwards) and outflow correctly computes the flow across it.</div><div><br></div><div>Yours,<br></div><div>Federico<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Il giorno lun 24 feb 2020 alle ore 16:20 Roland Haas <<a href="mailto:rhaas@illinois.edu">rhaas@illinois.edu</a>> ha scritto:<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 Federico,<br>
<br>
> I already track the SphericalSurface::sf_valid variable (line 570 of the<br>
> parfile). The sf_valid[2] ascii yields a -1 for each timestep, even after<br>
> t_merger = 1835.33. Yet, the other sf scalar outputs for the [2] surface<br>
> yield real numbers after t_merger (e.g., sf_max_radius[2] yields r =<br>
> 0.934217920021859 which is constant from the 13th iteration after the<br>
> merger onwards, and is *-nan* before t_merger).<br>
A value of "-1" means that (see the code snippet I had provided) that<br>
the horizon finder tried to look for a horizon and did not find it in<br>
that iteration.<br>
<br>
> As far as I can tell AHFinderDirect finds the 3rd horizon after the merger.<br>
> In my scalar output folder I have a <a href="http://BH_diagnostic.ah3.gp" rel="noreferrer" target="_blank">BH_diagnostic.ah3.gp</a> file which yields<br>
> real values starting from t_merger.<br>
If you see "sensible" values in the other fields but sf_valid is -1<br>
then this would mean that the horizon finder found the horizon at least<br>
once but is not currently finding it. Thus data on the surface if not<br>
guaranteed to actually be a good representation of the horizon since<br>
the "found it" could have been very long in the past.<br>
<br>
If you run eg<br>
<br>
git grep sf_max_radius<br>
<br>
in repos/einsteinanalysis/AHFinderDirect/src then you can see that<br>
sf_max_radius is only accessed in one place, namely in<br>
driver/BH_diagnostics.cc:<br>
<br>
--8<--<br>
681 // only try to copy AH info if we've found AHs at this time level<br>
682 if (! AH_data.search_flag) {<br>
683 sf_valid[surface_number] = 0;<br>
684 return;<br>
685 }<br>
686 <br>
687 // did we actually *find* this horizon?<br>
688 if (! AH_data.found_flag) {<br>
689 sf_valid[surface_number] = -1;<br>
690 return;<br>
691 }<br>
692 <br>
693 sf_active [surface_number] = 1;<br>
694 sf_valid [surface_number] = 1;<br>
695 sf_origin_x [surface_number] = this->origin_x;<br>
696 sf_origin_y [surface_number] = this->origin_y;<br>
697 sf_origin_z [surface_number] = this->origin_z;<br>
698 sf_mean_radius [surface_number] = mean_radius;<br>
699 sf_min_radius [surface_number] = min_radius;<br>
700 sf_max_radius [surface_number] = max_radius;<br>
--8<--<br>
<br>
and you can see that sf_max_radius is left at its old value if sf_valid<br>
is 0 (not trying to find) or -1 (not found).<br>
<br>
Thus if the horizon was ever found but is no longer then there will be<br>
sane values in sf_XXX but sf_valid will be -1.<br>
<br>
If you would like Outflow to continue using the (now invalid) horizon<br>
for a while until after it became invalid (for how long is up to you to<br>
decide) then you will have to add some grid scalar<br>
"horizon_last_valid_at" or so to Outflow (it must be checkpointed) and<br>
track when you last saw a "valid" horizon.<br>
<br>
Yours,<br>
Roland<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>