Dear all,
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 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.
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?
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.
Thank you very much for any help or insight you can provide!
Best,
Jordan Nicoules
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-217000do...) 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
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@mail.ubc.ca Sent: Wednesday, December 17, 2025 6:04:41 PM To: Jordan Nicoules Cc: users@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-217000do...) 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 .
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@mail.ubc.ca Sent: Wednesday, December 17, 2025 6:04:41 PM To: Jordan Nicoules Cc: users@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-217000do...) 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 .
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@mail.ubc.ca Sent: Tuesday, January 13, 2026 23:01 To: Jordan Nicoules Cc: users@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@mail.ubc.ca Sent: Wednesday, December 17, 2025 6:04:41 PM To: Jordan Nicoules Cc: users@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-217000do...) 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 .
Hello Jordan,
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):
Correct. The volume_forrm computed by Thornburg04 (and also Thornburg04nc I would guess) differs from the generic one computed otherswise.
This is something to be discussed with the Llama developers.
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) = detJinto
[... definition and attribution to h...] volume_form(i,j,k) = h(1)*h(2)*h(3) / abs(detJ)
Yes, a bug report would be most welcome.
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.
There is one remaining issue, namely that this will not take into account the case where there is overlap between coordinate patches. Thornburg04 handles this case by computing the overlap and reducing the volume_factor accordingly (0 for cells that are fully overlapped and a fraction of the volume if the overlap is not total).
I don't know on top of my head which other (functioning) coordinate system have overlap between patches.
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.
Well it would be step forward already to see the improvements that your changes bring.
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.
Well ideally the results of using Llama should be the same as not using Llama. Eg for a Cartesian only grid we currently have that the "sum" reduction in Carpet computes:
Sum(Grid_function) = integral(Grid_function)/(dx_coarse*dy_coarse*dz_coarse)
that is the "sum" reduction computes the Riemann integral of a grid function as if the coarsest grid spaccing was 1.0 instead of dx/dy/dz (since Carpet does not actually know that spacing).
So this should mean that for Llama the result for Cartesian coordinates should be the same and thus for non-Cartesian coordinates the result of the "sum" reduction should *also* be:
integral(Grid_function*volume_form)/(dx_coarse*dy_coarse*dz_coarse)
where now dx_coarse etc. are the grid spacing on some map, most conveniently a Cartesian map 0 if it exists (eg in Thornburg04 but not in Thornburg04nc).
Which is not quite what Llama produces right now, which would be:
Sum(Grid_function*volume_form) = int(Grid_function)
if I am not mistaken.
But that's a breaking change in Llama from its current behaviour, so will need some discussion.
Yours, Roland
Hi Roland,
Correct. The volume_forrm computed by Thornburg04 (and also Thornburg04nc I would guess) differs from the generic one computed otherswise.
This is something to be discussed with the Llama developers.
Thornburg04nc relies on the default behavior, which is why I find it a bit dangerous to have it inconsistent with Thornburg04. I think that Thornburg13 does it like Thornburg04. Maybe there are good reasons to have different behaviors for the volume form with different patch systems, but then again it forces outside code to check for the type of patch system.
Yes, a bug report would be most welcome.
Following your advice, I opened one earlier today.
I don't know on top of my head which other (functioning) coordinate system have overlap between patches.
Shouldn't that be done by the concerned patch system anyway, and currently give wrong results if it is not done?
But that's a breaking change in Llama from its current behaviour, so will need some discussion.
Sure, I understand. On that point, I already need to differentiate between Llama and non-Llama runs to fetch or not the volume form, so having different treatments was less of an issue.
Best,
Jordan
users@lists.einsteintoolkit.org
-
Jordan Nicoules -
Roland Haas