[Users] Safe check of variable Coordinates::volume_form_state

[Roland Haas]

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) = detJ
> 
> into
> 
>     [... 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

-- 
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 .