[Users] Safe check of variable Coordinates::volume_form_state
Jordan Nicoules
jnicoules at ua.pt
Fri Jan 16 11:11:46 CST 2026
Hi Roland,
Thank you very much for your answer and for pointing that out. Somehow I think that at the time, I was aware of this part of the source code, but I probably missed the fact that it is actually scheduled, even when the inverse Jacobian is not computed. So I guess that this really fixes my problem with checks and simplifies my life in that regard, because the parameter and flags are both set correctly then.
However, part of the reason why I thought this was not working properly, is because I was getting wrong numbers at the end of the day, but I attributed a wrong cause to it (bad memory access). Thanks to your comment about the spacing, and the piece of code you emphasized, I went back to doing some testing (with patch system Thornburg04nc), and I think that the problem lies in that the value assignment to the volume form is not consistent with what happens for Thornburg04 in the purely spherical part (with h the spacing):
// set volume form to deterimant of Jacobian
const CCTK_REAL det = fabs(( J11[ijk] * J22[ijk] * J33[ijk]
+ J12[ijk] * J23[ijk] * J31[ijk]
+ J13[ijk] * J21[ijk] * J32[ijk]
- J11[ijk] * J23[ijk] * J32[ijk]
- J12[ijk] * J21[ijk] * J33[ijk]
- J13[ijk] * J22[ijk] * J31[ijk]));
volume_form[ijk] = h[0]*h[1]*h[2]/det;
If it makes sense (and I'm not misled this time), I will instead open a bug report to change, in inverse-jacobian.F90,
volume_form(i,j,k) = detJ
into
[... definition and attribution to h...]
volume_form(i,j,k) = h(1)*h(2)*h(3) / abs(detJ)
which does indeed yield correct results for the volume integrals I'm computing. Also note the absolute value, because I've seen that detJ can be negative and it does impact the result at the end.
Again, I've limited my tests to Thornburg04 and Thornburg04nc patch systems (and only at t=0 for the latter), for the type of simulations I'm doing. So I don't know if it breaks things somewhere outside, and if the loss of backwards compatibility on that quantity may be a problem.
But I think this should improve consistency, and hopefully fit closer to the idea of the variable. This also prevents outsourcing that calculation with the spacing in any user thorn, hence limiting the information it needs to fetch from Coordinates.
Best,
Jordan
________________________________
From: Roland Haas <rhaas at mail.ubc.ca>
Sent: Tuesday, January 13, 2026 23:01
To: Jordan Nicoules
Cc: users at einsteintoolkit.org
Subject: Re: [Users] Safe check of variable Coordinates::volume_form_state
CUIDADO: Email de um sistema externo. Cuidado com links, anexos e pedidos de dados/senhas.
CAUTION: Email from an external system. Be careful with links, attachments, and requests for data/passwords.
Hello Jordan,
I took a quick look and Llama's schedule.ccl runs
Coordinates_SetInverseJacobian which sets volume_form_state and also
volume_form but only based on the Jacobian so does not take the grid
spacing into account or overlapping cells.
Thornburg04 then schedules its own routine afterwards which overwrites
those results.
volume_form_state = store_volume_form
[...]
if (store_volume_form /= 0) then
call calc_det3 (Jac, detJ)
! TODO: The volume form d^3x is defined as
! d^3x = da db dc * det(dx / da) * M,
! where M is the mask storing the nominal cell volume
volume_form(i,j,k) = detJ
end if
Yours,
Roland
> [CAUTION: Non-UBC Email]
>
> Hi Roland,
>
>
> Thank you very much for your message! I had seen the shorter version
> of your reply in the minutes indeed, but thank you for giving a more
> detailed answer!
>
>
> When I come back from the holiday season, I will gladly file the bug
> report. I would also say that this should be straightforward to fix,
> at least about initializing the variable to 0 by default.
>
>
> The consistency between setting the parameter
> Coordinates::store_volume_form=yes and having the code still set
> Coordinates::volume_form_state = 0 depending on the actual
> implementation of the volume form -- if it's even needed -- might be
> a different question though. For instance, I think that technically
> the system Thornburg04nc should have the volume form computed, but it
> doesn't. I suppose that the consistency could be handled at
> ParamCheck for example. Please tell me if that would require another
> subsequent bug report or not.
>
> I would say that this renders the usage of CCTK_VarDataPtr on
> Coordinates::volume_form still insufficient because its storage
> relies on the parameter value, although I get your point about
> abusing poisoning there. However (from the top of my head), I
> remember that when I used poisoning with Llama (and MLBSSN), the
> norm2 of the Hamiltonian constraint was NaN, so I had to turn it off.
>
>
> Let me know if you'd like me to help on the implementation. I'm not
> sure about how the process goes and what is the corresponding
> involvement, so I don't want to make promises, but I'd happy to
> contribute to these (hopefully easy) fixes, if it saves you more time
> than it wastes.
>
>
> Best wishes for the holiday season!
>
>
> Jordan
>
>
>
> ________________________________
> From: Roland Haas <rhaas at mail.ubc.ca>
> Sent: Wednesday, December 17, 2025 6:04:41 PM
> To: Jordan Nicoules
> Cc: users at einsteintoolkit.org
> Subject: Re: [Users] Safe check of variable
> Coordinates::volume_form_state
>
>
> CUIDADO: Email de um sistema externo. Cuidado com links, anexos e
> pedidos de dados/senhas. CAUTION: Email from an external system. Be
> careful with links, attachments, and requests for data/passwords.
>
> Hello Jordan,
>
> Sorry for the delayed response. This dropped off my radar after
> agreeing to try and respond during the last ET call.
>
> > I'm trying to extend our analysis thorn to tackle multipatch grids
> > built with Llama. I want to use Coordinates::volume_form in the
> > computation of volume integrals, by setting the corresponding
> > parameter Coordinates::store_volume_form = yes.
> >
> >
> > Our standard use case relies on Coordinates::coordinate_system =
> > "Thornburg04", for which it works like a charm. However, except for
> > this system and the analogous Thornburg13 (I have not tested it so
> > far), it seems the variable Coordinates::volume_form_state is not
> > modified nor initialized, whatever the value of the parameter
> > Coordinates::store_volume_form.
>
> That would be a bug in the Llama coordinate systems. They should all
> at least set up the variable correctly. Would you mind filing a bug
> report using
>
> https://bitbucket.org/einsteintoolkit/tickets/issues/new
>
> > That makes the parameter insufficient
> > by itself to perform checks. Moreover, since the flag is not
> > initialized, we also can't simply check its value, which can be
> > polluted. I have tried to use Carpet::poison_new_timelevels = yes
> > and CarpetLib::poison_new_memory = yes, but it didn't seem to
> > solve the issue.
>
> You can check if a variable has storage using the CCTK_VarDataPtr
> function call
> (https://einsteintoolkit.org/referencemanual/ReferenceManual.html#x1-217000doc)
> but it won't tell you if the values found in there are sensible
> (admittedly you could abuse poisoning or this).
>
> > Hence, I would like to ask if there's another way to check for
> > uninitialized values, a safe way to use the variables I mentioned in
> > a general case, or any information/subtlety I may have missed?
>
> I think you captured the intended behaviour and the bugs present.
> Should be (one hopes) easy to fix at least in that all coordinate
> systems should indicate correctly if the volume form is set up or now.
>
> > Ultimately, I can directly check Coordinates::coordinate_system, and
> > output a warning/error if it's not of the types mentioned above.
> > That feels a bit unsatisfactory, although I may default to it.
>
> Yes, looking at another thorns parameter is somewhat error prone, in
> particular if it is a parameter that is private and thus could change
> meaning at any moment.
>
> Yours,
> Roland
>
> --
> My email is as private as my paper mail. I therefore support
> encrypting and signing email messages. Get my PGP key from
> http://pgp.mit.edu .
--
My email is as private as my paper mail. I therefore support encrypting
and signing email messages. Get my PGP key from http://pgp.mit.edu .
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.einsteintoolkit.org/pipermail/users/attachments/20260116/6c503629/attachment.htm>
More information about the Users
mailing list