Conversation
AsterX: add code for RePrimAnd C2P
| pv.vel = (mu * cache.x) * (r_u + (rb * mu) * b_u); | ||
| pv.w_lor = cache.w; | ||
| } | ||
|
|
There was a problem hiding this comment.
Exact same comment as for Palenzuela, see here:
https://github.com/EinsteinToolkit/AsterX/pull/126/changes#r2886212547
MChabanov
left a comment
MChabanov
left a comment
There was a problem hiding this comment.
Some comments from my side! Code looks very good!
| // where r^i = S^i / D, b^i = B^i / sqrt(D), (r.b) = (B^i S_i)/(D^(3/2)). | ||
| // ------------------------------------------------------------------ | ||
|
|
||
| const vec<CCTK_REAL, 3> mom_up = calc_contraction(gup, cv.mom); |
There was a problem hiding this comment.
Just an idea: One can extract the code below, where prims are computed, to a separate function, similar to our other C2Ps. Eg for Palenzuela see xPalenzuelaToPrim in
diagnostic, add vwlim condition in Palenzuela C2P
MChabanov
left a comment
MChabanov
left a comment
There was a problem hiding this comment.
Great effort! I'm very impressed. Most comments are minor and are mainly related to consistency with the old code. Thanks!
| (void)status; | ||
| return; | ||
| } | ||
| rep.set_soft_root_conv(); |
There was a problem hiding this comment.
line 409: the root finding finds the root for rho, such that max(1, max(rho_l, rho_r) ) might always evaluate to 1 ... -> a better choice for the lower limit of the scale might be just 0, ie
Replace with:
fmax(CCTK_REAL(0.0), fmax(abs(result.first), abs(result.second)));
| return; | ||
| } | ||
| rep.set_soft_root_conv(); | ||
| } |
There was a problem hiding this comment.
line 494: similar comment as for the entropy C2P.
Palenzuela root finds for hW which should always be >= 1.0. So, that's fine. However, some EOS tables have negative eps in their tables which might lead to a h < 1 in the code. This means that formally you would want h_min (from the table) here as the lower bound and having h>=1 would mess with these values ...
| pv.temperature = eos_3p->temp_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); | ||
| pv.entropy = eos_3p->kappa_from_rho_eps_ye(pv.rho, pv.eps, pv.Ye); | ||
| rep.adjust_cons = true; | ||
| } |
There was a problem hiding this comment.
Velocity and temperature are limited by the code implemented in the function
c2p::prims_floors_and_ceilings
called below in line 537. Additionally, limiting eps or temp should be equivalent in AsterX to my understanding, hence we might already get away without these limits here.
I suggest to explore whether the limiting we already have is sufficient and if not, this code could be moved to the function c2p::prims_floors_and_ceilings in a consistent way. To avoid calling the EOS too often and having code duplication.
| atmo.set(pv, cv, glo); | ||
| return; | ||
| } | ||
|
|
There was a problem hiding this comment.
Note that this functionality already exists in line 226 of con2prim.cxx
AsterX/AsterX/src/con2prim.cxx
Line 226 in 19a463e
| cv.DYe); | ||
| set_to_nan(pv, cv); | ||
| return; | ||
| } |
There was a problem hiding this comment.
isfinite(cv.dens) is checked below in line 116. cv.dens <= 0 should be captured before
AsterX/AsterX/src/con2prim.cxx
Line 226 in 19a463e
Hence, I'm not sure if we need this check here?
| CCTK_REAL Ye_raw = cv.DYe / cv.dens; | ||
| const CCTK_REAL Ye = fmin(fmax(eos_3p->rgye.min, Ye_raw), eos_3p->rgye.max); | ||
| const bool ye_clipped = | ||
| (Ye_raw < eos_3p->rgye.min) || (Ye_raw > eos_3p->rgye.max); |
There was a problem hiding this comment.
That's not part of this PR but this line makes me think that we might want to do something similar for Palenzuela ...
| if (std::isfinite(a) && std::isfinite(b) && b >= a) { | ||
| const CCTK_REAL width = b - a; | ||
| const CCTK_REAL scale = | ||
| fmax(CCTK_REAL(1.0), fmax(std::abs(a), std::abs(b))); |
There was a problem hiding this comment.
Same comment as before for Palenzuela and Entropy:
To my understanding RPA solves for 1/(hW) which has the limits 0 < mu < 1/hmin ...
That means you will likely choose 1 for the scale ...
I would be good to double-check that, in case that's true the better lower limit would be again 0.0, so, replace line 181 with
fmax(CCTK_REAL(0.0), fmax(std::abs(a), std::abs(b)));
| { | ||
| const auto rgeps = eos_3p->range_eps_from_rho_ye(pv.rho, pv.Ye); | ||
| if (cache.eps_raw < rgeps.min || cache.eps_raw > rgeps.max) { | ||
| rep.adjust_cons = true; |
There was a problem hiding this comment.
Why are you setting the flag for recomputation if eps didn't change here?
2 participants
Add RPA support to Con2PrimFactory.
Code tested on Vista: