improve diffusion-dislocation creep - #7270
Conversation
0286ae9 to
96e742e
Compare
alarshi
left a comment
alarshi
left a comment
There was a problem hiding this comment.
Thank you for making this PR! I have wondered in the past why we allow n != 1 diffusion creep stress exponents. I just have a small clarifying comment, otherwise it looks good to me. I don't know whether we use a value other than 1 for the diffusion creep stress exponents in any of the tests, but we will find out.
| // Our starting guess assumes that all strain is accommodated by diffusion creep. | ||
| // If the diffusion creep prefactor is very small, that means that the diffusion viscosity is very large. | ||
| // In this case, use the maximum viscosity instead to compute the starting guess. | ||
| double stress_ii = (prefactor_stress_diffusion > (0.5 / maximum_viscosity) |
There was a problem hiding this comment.
Shouldn't this line also be modified to stress_ii = (prefactor_stress_diffusion > 2 * maximum_viscosity * edot_ii ?
There was a problem hiding this comment.
No, but your comment made me realise that this logic was very unclear, so I have refactored it. prefactor_stress_diffusion was equal to edot_ii/stress_ii. I think the new code is much clearer and the impact on performance will be negligible.
| // Power law creep equation (grain size independent): | ||
| // edot_ii = A * stress^n * exp(-(E + P*V)/(RT)) | ||
| // By taking the logarithm: | ||
| // log(edot_ii) = std::log(A) + n*std::log(stress) - (E + P*V)/(RT) |
There was a problem hiding this comment.
I really liked reading the explanation, thank you for expanding the comments so clearly.
Could we maybe modify this for consistency?
log(edot_ii) = log(A) + n*log(stress) - (E + P*V)/(RT) or std::log(edot_ii) = std::log(A) + n*std::log(stress) - (E + P*V)/(RT) ?
There was a problem hiding this comment.
Good idea. I've now modified this.
791d797 to
d7dbefb
Compare
|
@alarshi good for another look :) |
|
/rebuild |
gassmoeller
left a comment
gassmoeller
left a comment
There was a problem hiding this comment.
Thanks, this makes a lot of sense. Please check my two comments. I would phrase the documentation of the parameter stricter, there is no reason to keep around unused parameters in the long run.
If you implement my suggestions, this should be ready.
| "The stress exponent for diffusion creep is almost always equal to one. " | ||
| "If only one value is given, then all use the same value. Units: None."); | ||
| "If only one value is given, then all use the same value. " | ||
| "This is a redundant parameter, as the stress exponent for diffusion creep must always equal one. " |
There was a problem hiding this comment.
| "This is a redundant parameter, as the stress exponent for diffusion creep must always equal one. " | |
| "This is a deprecated parameter, as the stress exponent for diffusion creep must always equal one. This parameter will be removed in a future version." |
| Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 2.19823e-05 | ||
| Relative nonlinear residual (total system) after nonlinear iteration 2: 2.19823e-05 | ||
| Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.68721e-16, 1.05609e-16, 2.19834e-05 | ||
| Relative nonlinear residual (total system) after nonlinear iteration 2: 2.19834e-05 |
There was a problem hiding this comment.
It doesnt look like these two tests are outputting much diagnostic output. This makes it hard to judge if this change happened because of a rounding difference, or if something substantial changed. Could you add either the material statistics postprocessor output, or the depth average postprocessor output to the output folders of these tests? Just so we can check how large the difference is in these outputs?
d7dbefb to
d5a715f
Compare
78592ef to
9020be0
Compare
aeec0cc to
356a628
Compare
|
#7279 has been merged, so we can now directly compare main to this PR. The test results have changed a bit. I believe that this is due to the changed starting guess for the stress for the Newton solve: The only other bit of code that has changed is the calculation of the diffusion viscosity, but the actual mathematics are the same, I just stripped out terms involving the stress exponent: https://github.com/geodynamics/aspect/pull/7270/changes#diff-c305db66b306a8a86c9dd166a2a8eb4c0d3fe56681147d2266c672cc6148dd9eR142-R150 So I think this is ready to merge. |
3 participants
This PR:
diffusion_creep.cc,dislocation_creep.ccanddiffusion_dislocation.ccstress_exponentsfrom processed parameters indiffusion_creep.cc- these were never used in the rheological model, and only partially used indiffusion_dislocation.cc, making the (non)implementation unclear.diffusion_dislocation.cc. This should have negligible effect on performance or accuracy, but it makes sense to use the correct value.