model hangs or floating point exception error in visco-plastic module

mibillen commented 3 years ago

I came across an error while creating my visco-elasto-plastic subduction models. In optimized mode the model just hangs. In debug mode I get some kind of floating point exception error (SIGFPE) sometimes with no information about where this error is occurring. In changing to debug mode and trying to reduce the convergence paramters to reproduce the error more quickly, I then got this error message: *** Timestep 3: t=383.779 years, dt=160.397 years Solving mesh velocity system... 18 iterations.

TimerOutput objects finalize timed values printed to the screen by communicating over MPI in their destructors. Since an exception is currently uncaught, this synchronization (and subsequent output) will be skipped to avoid a possible deadlock.

Exception 'ExcMessage ("All parameter values must be greater than 0 for harmonic averaging!")' on rank 0 on processing:

An error occurred in line <813> of file </Users/billen/Box-Sync/Aspect-Src/aspect/source/material_model/utilities.cc> in function double aspect::MaterialModel::MaterialUtilities::average_value(const std::vector &, const std::vector &, const enum CompositionalAveragingOperation &) The violated condition was: parameter_values[i] > 0 Additional information: All parameter values must be greater than 0 for harmonic averaging!

Aborting!

This error has to do with the harmonic average (compositional_averaging) which checks that the viscosity is greater than 0 (material_model/utilities.cc average_value function on about line 810 - here I found that there was a viscosity value was nan). When I did I managed to track down why the viscosity was NaN (using "print" statements in visco_plastic.cc), I found the error originated to line 179 of .../rheology/visco_plastic.cc where the composite viscosity is calculated from the diffusion creep and dislocation creep. The problem was that I had a temperature value of only 88 K (due to too-high velocity with the free surface I think) and the diffusion creep viscosity was about 1.5e104 while the dislocation creep was about 1e210. So, the top line of the composite_viscosity (eta_df*eta_ds) = 15e314 ended up being larger than a double. It might be useful to have a check here on log10(eta_df)+log10(eta_ds) > 300. Could also be useful to have an assert before calculating the last line of the (bounding the viscosity between visc_min and visc_max), to make sure the viscosity is not a NaN before this last step.

bangerth commented 3 years ago

Good detective work! I think adding such an assertion would make good sense. Would you like to give such a patch a try?

bangerth commented 3 years ago

Separately, I think it would be useful to compare against something like std::log10(std::numeric_limits<double>::max() / 2 instead of 300. This makes it clear to the reader that you're not trying to check against an arbitrary threshold, but that you're concerned about overflow.

mibillen commented 3 years ago

Thanks Wolfgang, I figured there was a better way to check this, but I would not have known what it is. Incidentally, I kind of new to look for this because I needed the same kind of check in CItcom … back in the day ;-), but that was when we were just using floats, so the low temperatures in the viscosity laws (even when not an overshoot) would exceed the largest float.

I’ll give the patch a try - but I need to submit my proposal first, so this will need to wait a week or so. Magali

On Jan 25, 2021, at 5:30 PM, Wolfgang Bangerth notifications@github.com wrote:

Separately, I think it would be useful to compare against something like std::log10(std::numeric_limits::max() / 2 instead of 300. This makes it clear to the reader that you're not trying to check against an arbitrary threshold, but that you're concerned about overflow.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3959#issuecomment-767222881, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACS64UWOBZBCAWTFVGDJYDDS3YLL5ANCNFSM4WSNG2ZQ.

Professor of Geophysics Chair, Geology Graduate Program Geophysics Minor Advisor

Department of Earth and Planetary Sciences Univ. of California, Davis Davis, CA 95616 2129 Earth & Physical Sciences Bldg. Office Phone: (530) 752-4169 https://magalibillen.faculty.ucdavis.edu/

Pronouns: She/her/hers First Generation College Student: https://firstgen.ucdavis.edu/

elodie-kendall commented 3 years ago

Thanks a lot for this Magali! I look forward to hearing if your fix works. Elodie

mibillen commented 3 years ago

Hi Becky and Trey,

We are going to be using the terminal a lot on peloton. So here is a 50 minute intro tutorial to the linux/unix terminal, https://ubuntu.com/tutorials/command-line-for-beginners#3-opening-a-terminal https://ubuntu.com/tutorials/command-line-for-beginners#3-opening-a-terminal

If you’ve used a terminal before it might just be a refresher with some handy hints. Otherwise, it starts at the beginning so you can have a clear understanding of how to navigate, access files, etc… Peloton uses ubuntu, but the mac terminal is 95% the same.

Magali

Professor of Geophysics Chair, Geology Graduate Program Geophysics Minor Advisor

Department of Earth and Planetary Sciences Univ. of California, Davis Davis, CA 95616 2129 Earth & Physical Sciences Bldg. Office Phone: (530) 752-4169 https://magalibillen.faculty.ucdavis.edu/

Pronouns: She/her/hers First Generation College Student: https://firstgen.ucdavis.edu/

elodie-kendall commented 3 years ago

Any luck Magali? I'm trying myself but I'm new to ASPECT so it's taking a bit of time. Thanks :)

mibillen commented 3 years ago

I’m sorry, but I wont’ have time to try this for several weeks due to administration and limited work time wtih children home during the pandemic.

On Feb 4, 2021, at 5:03 AM, elodie-kendall notifications@github.com wrote:

Any luck Magali? I'm trying myself but I'm new to ASPECT so it's taking a bit of time. Thanks :)

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/geodynamics/aspect/issues/3959#issuecomment-773288287, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACS64USLCVVPMLSDY2652PLS5KLIZANCNFSM4WSNG2ZQ.

Professor of Geophysics Chair, Geology Graduate Program Geophysics Minor Advisor

Department of Earth and Planetary Sciences Univ. of California, Davis Davis, CA 95616 2129 Earth & Physical Sciences Bldg. Office Phone: (530) 752-4169 https://magalibillen.faculty.ucdavis.edu/

Pronouns: She/her/hers First Generation College Student: https://firstgen.ucdavis.edu/

geodynamics / aspect

model hangs or floating point exception error in visco-plastic module #3959

TimerOutput objects finalize timed values printed to the screen by communicating over MPI in their destructors. Since an exception is currently uncaught, this synchronization (and subsequent output) will be skipped to avoid a possible deadlock.

Aborting!