boutproject
Laplace optimisations when metric coefficients are constant
Running e.g. examples/bout_runners_example/3D_diffusion:
'Laplace( Field3D ) Sorry, D2DXDY not implemented'
This error appeared after git pull this morning (do not know tag of latest working version, sorry)
This is because D2DXDY used to just return zero, but now throws this error message. The possibilities are:
a) The metric tensor components g12 (x-y) are zero for all coordinates used in practice (as far as I know), so we could remove this term from the Laplace operator. This might then surprise someone later if they want g12 != 0
b) We could implement D2DXDY, which would currently require taking derivatives in one direction, communicating, and taking derivatives the other way. Most of the time this work would be done, and then set to zero anyway.
Would it be useful to give metric tensors a flag (in Field2D?) to specify that they are zero. This could then be used in operators like Laplace to decide whether to evaluate a term.
bool Field2D.isNonZero(); Return true if not zero
then have functions in Field2D
void Field2D.markZero(); // Mark this field as zero. The code could check values, and throw an error if not zero, or maybe set itself to zero.
void Field2D.markNonZero(); // Set back to non zero. This would be the default
Laplace should do exactly a Laplace operation so solution a) opens up for future bugs. Thererefore, I vote for a solution in line with Ben solution "c"):
Would it be useful to give metric tensors a flag (in Field2D?) to specify that they are zero. This could then be used in operators like Laplace to decide whether to evaluate a term.
bool Field2D.isNonZero(); Return true if not zero
then have functions in Field2D
void Field2D.markZero(); // Mark this field as zero. The code could check values, and throw an error if not zero, or maybe set itself to zero.
void Field2D.markNonZero(); // Set back to non zero. This would be the default
This kind of solution could also save a small amount of computational time in shifted metric cases where g^{xz} is also zero but is calculated in e.g. 'Laplace()'
Concerning a) It would be possible to check if g12 is zero, and throw an error if not ... it still needs implementing in the future ....
Concerning b)
It would be possible to communicate the corner cells without much overhead. At the X-point there is more than one option what to communicate, but then the D2DXDY has also this issue, as near the X-Point the symmetry of second derivatives is broken ...
Concerning markZero: Instead of having just the option to set a field to zero or not - it might be beneficial to be able to set the whole field to a constant value. I don't think it is significantly more complicated than setting to zero or not, as that has still the trouble if I want to set an value of an field which should be constant - should it set the whole field to this value? Should it throw an error? Should it silently convert to non-constant? I guess throwing is the best option. The benefit would be, I could also set to one - and multiplications would be much faster ... Also multiplying two constant fields would just need a single multiplication ...
I have thought about trying to give fields a single constant value, but think it would be more complicated to implement. Code like Laplace() would have to check if it was constant, then if that constant was zero. I agree though that it could be a nice optimisation in some cases.
A field could have functions like
void Field2D.setConstant(BoutReal value); bool Field2D.isConstant(); // Return true if a constant value BoutReal Field2D.getConstant(); // Get the constant value
I think errors could be minimised for this solution or the setting to zero flag by:
- When setConstant is called (or markZero), all the field values are set to that constant. If the field is accessed using existing code, it will see the correct value.
- If any value is assigned to the field (operator= called), the constant flag becomes false. Some optimisation would be possible, like operator=(BoutReal) could remain constant, and assignment from a constant field could also be constant.
Perhaps we could have the above functions, then also
bool Field2D.isNotZero();
internally, the flag this returns could be decided in setConstant() e.g.
class Field2D { private: bool is_constant; // Is this field a constant value? bool is_nonzero; // Is this field nonzero?
public: void setConstant(BoutReal val) { is_constant = true; if(fabs(val) < 1e-15) { is_nonzero = false; }else { is_nonzero = true; } // Set all internal data to the constant value data = val; } }
The important part is to make sure that these flags degrade as soon as the data is modified in any way. This is easy for assignment operations, but difficult for things like
Field2D f;
f.setConstant(0.0); // Set to a constant
f(2,3) = 1.0; .// Modify an element // Error: f still constant!
Unfortunately I think fixing this would mean putting code into operator(), which would be on the inner loop of iterations. Setting to zero (or a constant) would have to be used with care, or it could lead to unexpected results.
As we return a reference to a value, we do not know whether this value is written or read. Otherwise we could just raise an error if the value is written (in case CHECK is enabled)
On the other hand we allow the user to break stuff. This is nothing new, and I don't think we should prevent this, only make it not to easy to do it. With having the .check() function, we could amend it to also check the values, if it should be constant.
Thereby this flag would be an promise to the library to keep this field constant. If this promise is broken, unexpected behavior is to be expected.
I think PR #263 fixed this issue, so perhaps the optimisation discussion is an enhancement request