ianhinder / Kranc

A Mathematica package for generating code for solving time dependent partial differential equations
http://kranccode.org
GNU General Public License v2.0
29 stars 10 forks source link

Suspicious Jacobian code #106

Closed eschnett closed 10 years ago

eschnett commented 10 years ago

This generated code looks suspicious: {{{ const bool use_jacobian1 = (!CCTK_IsFunctionAliased("MultiPatch_GetMap") || MultiPatch_GetMap(cctkGH) != jacobian_identity_map) && strlen(jacobian_group) > 0; }}}

It says that, if MultiPatch_GetMap is not aliased, then we assume we should use a Jacobian? The opposite should be true, since we then won't be using a multi-block system.

barrywardell commented 10 years ago

I recently noticed this line too and was confused by it. I concluded it was probably fine as long as you don't try to use single-block and set the jacobian_group parameter.

It's saying that we should use a Jacobian if jacobian_group is defined and either we're using single-block (MultiPatch_GetMap is not aliased) or we're using multi-block with a Jacobian which is not the identity. If it never makes sense to use a Jacobian with single-block, then I agree this is suspicious and should probably change to something like (is the strlen piece needed then?):

(CCTK_IsFunctionAliased("MultiPatch_GetMap") && MultiPatch_GetMap(cctkGH) != jacobian_identity_map)
&& strlen(jacobian_group) > 0
ianhinder commented 10 years ago

This is confusing, yes, but I think it is correct. GenericFD has a jacobian defined in its interface.ccl so that tests can have a jacobian without having to have Llama (which was a private code at the time, and is still not a part of the Einstein Toolkit). What this is saying is that if there is no Llama available but there is a Jacobian group specified in the parameter file, then we should use it. This Jacobian used for testing has no storage assigned by default. It could be renamed as "test_jacobian" or something to make this a bit clearer. It should also be moved into the generated thorn.

eschnett commented 10 years ago

The code seems correct after all.