Do not solve advection system if RHS is zero #898
Conversation
|
👍 I just realized we can use |
| build_advection_preconditioner(AdvectionField::temperature(), | ||
| T_preconditioner); | ||
| solve_advection(AdvectionField::temperature()); | ||
| if (system_rhs.block(AdvectionField::temperature().block_index(introspection)).l2_norm() <= std::numeric_limits<double>::min()) |
There was a problem hiding this comment.
This should just compare to numeric zero, for two reasons.
First, x <= numeric_limits::min() is true if either x==numeric_limits::min() (which is a floating point comparison for equality that is unlikely to ever be true), or if x==0. The definition of numeric_limit::min() is so that the condition x < numeric_limit::min() is never true.
Second, rhs.l2_norm() has physical units, whereas numeric_limits::min() does not, and so a comparison cannot be correct without a scaling factor.
There was a problem hiding this comment.
This is not quite correct (there is -0.0 and denormalized values that I catch this way), see http://cpp.sh/2nbma for example.
Second, rhs.l2_norm() has physical units,
Yes, we should probably do something like 1e-16*something_sensible_with_sensible_units, but I am not confident how to do this robustly so we went for the easy option.
There was a problem hiding this comment.
So in which cases do you think you could get -0.0 or a denorm out of the l2 norm? I think that's an academic consideration. Just compare against zero and be done with it. That's exactly the case you wanted to catch after all.
There was a problem hiding this comment.
I am not sure it is needed, but with your suggestion I get
a.cc: In function ‘int main()’:
a.cc:14:26: warning: comparing floating point with == or != is unsafe [-Wfloat-equal]
std::cout << "" << (a==0.0);
with g++ 6.1. Any suggestions?
There was a problem hiding this comment.
Baa, what a stupid warning! Comparing against 0 should definitely be safe.
I guess in that case leave it as it is, and add a comment to the effect that we really just want to compare against zero but can't because of the warning...
There was a problem hiding this comment.
Okay, I added the comment, and after I noticed that I would have to add it in about 20 places, I moved all my new additions into solve_advection so that we do not have to duplicate the code.
|
Nice improvement. There are lots of failing tests right now. I can't seem to connect to @tjhei 's tester at the moment, so can't say what's going on, but I suspect that for lots of tests the check triggers. Can you update the test outputs? |
|
@tjhei -- what are the circumstances where the matrix is zero? This would make no sense in a time stepping scheme... |
|
Yes, exactly, but right now if the matrix is zero and the RHS is nonzero, we trigger some error in the advection preconditioner, and it would take the user some time to figure out what's going on. So we thought, maybe we could just add an |
|
Yes, I agree that the assertion is a good idea. I just can't come up with any scenario in which the matrix would be zero. (Though there are now many more possibilities to screw up with the new assembler structure.) |
|
So we thought if c_p, thermal conductivity, and the entropy viscosity stabilization are zero, the matrix will be zero. (I am not really sure in which cases exactly the entropy viscosity stabilization is zero, maybe when the field itself is zero, but I had a model where the matrix was the problem, and I got a not very useful exception in the advection preconditioner.) There are terms on the RHS that are independent of that (e.g. the radiogenic heating), so we could end up with a zero matrix and a nonzero RHS. |
Clemson had a planned power outage this weekend (yes, seriously). Should be back to normal later today. |
jdannberg
force-pushed
the
zero_RHS
branch
3 times, most recently
from
1d70fc2
to
1104082
Compare
|
I see. Yes, if |
jdannberg
force-pushed
the
zero_RHS
branch
3 times, most recently
from
44e7057
to
884f8ce
Compare
|
This looks good to me. @tjhei -- there is a failure with this test here: http://www.math.clemson.edu/~heister/cib/aspect-8.4.1/cc22a221934bec4c0d561d798af347162a9d0c94/changes-clang.diff |
I was going to fix this together with Juliane later today. I think it is a very bad idea merging pullrequests that do not pass all tests. |
|
Yes, my bad. I had looked at the error message, found it unrelated to the current patch, and decided that we can handle these errors in patch review. I knew you were eventually going to fix this, but you're right that I could also just have waited for that to happen. |
No problem. I guess as long as we try to keep up with failing tests, we will be fine. ;-) |
Right now if the matrix and the right-hand side are zero (for example for the temperature advection), we get a floating point exception in the advection preconditioner that is not really helpful to figure out what the problem is. In addition, if there is a problem that does not need the temperature, but needs compositional fields, or it has a non-zero temperature/composition only from a given time step on, then we do not really want to solve the advection system in the time steps before.
This pull request changes the solver schemes, so that in case of a zero RHS we do not solve for the advection field, but just set it to zero everywhere.