[aspect-devel] muparser and Floating point exception from divide by zero.
Bob Myhill
myhill.bob at gmail.com
Fri Jun 29 12:16:36 PDT 2018
Dear Magali,
You're using a muParser if syntax that was deprecated a few years ago. The
new syntax uses lazy evaluation to evaluate only the required branch of the
if-statement. Try this:
x>0 ? TS+(TM-TS)*(1-erfc((YMAX-y)/(2*sqrt(K*x/VSUB)))) : TM
You might also want to consider an eps in your if-condition, otherwise very
small numbers might also give you floating point exceptions.
Let me know if this doesn't work.
Best wishes,
Bob
On 29 June 2018 at 09:43, Magali Billen <mibillen at ucdavis.edu> wrote:
> Hello all,
> I have small mystery I could use some help with.
>
> I’m using a function to define the initial temperature field which is only
> defined for x (horizontal position)
> greater than zero because x shows up in the denominator. At x=0 the
> temperature equals a constant value, TM
>
> If I add a variable EPS=1e-12, and use (x+EPS), this fixes the problem.
> However, if I instead use an if statement: if (x>0,..., TM), then I get a
> floating point exception, as I did in the case with no (x+EPS).
>
> Any ideas why the if statement doesn’t work? This suggest to me that the
> equation is actually evaluated everywhere and then the if
> condition is used, rather than the if condition being checked first (as it
> would be if written in a c program). Perhaps this is how
> the muparser works, but this seems quite strange?
>
> I’d like to use the if statement, because for my actual runs, I have a
> complicate temperature function, and it would be cleaner (easier to debug),
> not to have a bunch of (x+EPS) in the equations.
>
> I’ve pasted in the two different functions calls below. The errors are
> pasted in below that.
> And, I’ve attached the prm-file (with the if statement).
>
> Here are the actual function calls:
>
> ==> Option 1 - with eps = 1e-12 (this works, no floating point exception)
>
> subsection Initial temperature model
> set Model name = function
>
> subsection Function
> set Variable names = x,y
> set Function constants = EPS=1e-12, YMAX=1e6, TS=0.000e+00,
> TM=1.400e+03, \
> K=1.000e-06, VSUB=1.585e-09
> set Function expression = TS+(TM-TS)*(1-erfc((YMAX-y)/(
> 2*sqrt(K*(x+EPS)/VSUB))))
> end
> end
>
> ==> Option 2 with if statement. (this does not work, I get a floating
> point exception)
>
> subsection Initial temperature model
> set Model name = function
>
> subsection Function
> set Variable names = x,y
> set Function constants = YMAX=1e6, TS=0.000e+00, TM=1.400e+03, \
> K=1.000e-06, VSUB=1.585e-09
> set Function expression = if(x>0, TS+(TM-TS)*(1-erfc((YMAX-y)/(
> 2*sqrt(K*x/VSUB)))),TM)
> end
> end
>
>
> > more *.output
> ------------------------------------------------------------
> -----------------
> -- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
> -- . version 2.1.0-pre (master, f5e6a41d5)
> -- . using deal.II 9.0.0
> -- . using Trilinos 12.10.1
> -- . using p4est 2.0.0
> -- . running in DEBUG mode
> -- . running with 8 MPI processes
> -- How to cite ASPECT: https://aspect.geodynamics.org/cite.html
> ------------------------------------------------------------
> -----------------
>
> Number of active cells: 1,280 (on 5 levels)
> Number of degrees of freedom: 17,316 (10,626+1,377+5,313)
>
> [c12-7:17151] *** Process received signal ***
> SIGFPE received
> [c12-7:17151] Signal: Floating point exception (8)
> [c12-7:17151] Signal code: Floating point divide-by-zero (3)
> [c12-7:17151] Failing at address: 0x7f37f5fc163a
> [c12-7:17151] [ 0] /lib/x86_64-linux-gnu/libpthread.so.0(+0x12890)[
> 0x7f37fac91890]
> [c12-7:17151] [ 1] /usr/lib/x86_64-linux-gnu/libmuparser.so.2(_
> ZNK2mu10ParserBase16ParseCmdCodeBulkEii+0x24a)[0x7f37f5fc163a]
> [c12-7:17151] [ 2] /usr/lib/x86_64-linux-gnu/libmuparser.so.2(_
> ZNK2mu10ParserBase11ParseStringEv+0x29)[0x7f37f5fc8999]
> [c12-7:17151] [ 3] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZNK6dealii14FunctionParserILi2EE5valueERKNS_5PointILi2EdEEj+0x2c
> 2)[0x7f3803f1864c]
> [c12-7:17151] [ 4] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZNK6dealii9Functions14ParsedFunctionILi2EE5valueERKNS_5PointILi2
> EdEEj+0x10)[0x7f3803f85060]
> [c12-7:17151] [ 5] /home/billen/AspectProjects/aspect/build/aspect(_
> ZNK6aspect18InitialTemperature8FunctionILi2EE19initial_
> temperatureERKN6dealii5PointILi2EdEE+0x6d)[0x5
> 5b227f7fe15]
> [c12-7:17151] [ 6] /home/billen/AspectProjects/aspect/build/aspect(_
> ZNK6aspect18InitialTemperature7ManagerILi2EE19initial_
> temperatureERKN6dealii5PointILi2EdEE+0x47)[0x55
> b227f83f33]
> [c12-7:17151] [ 7] /home/billen/AspectProjects/
> aspect/build/aspect(_ZNSt17_Function_handlerIFdRKN6dealii5PointILi2
> EdEEESt5_BindIFMN6aspect18InitialTemperature7ManagerILi
> 2EEEKFdS4_ESt17reference_wrapperISA_ESt12_PlaceholderILi1EEEEE9_M_
> invokeERKSt9_Any_dataS4_+0x20)[0x55b228455ceb]
> [c12-7:17151] [ 8] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZNKSt8functionIFdRKN6dealii5PointILi2EdEEEEclES4_+0xe)[0x7f3803e
> b3b52]
> [c12-7:17151] [ 9] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZNK6dealii38VectorFunctionFromScalarFunctionObjectILi2EdE12vecto
> r_valueERKNS_5PointILi2EdEERNS_6VectorIdEE+0x4a)[0x7f3803ecee64]
> [c12-7:17151] [10] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZNK6dealii8FunctionILi2EdE17vector_value_listERKSt6vectorINS_5Po
> intILi2EdEESaIS4_EERS2_INS_6VectorIdEESaISA_EE+0xa4)[0x7f3803eb835c]
> [c12-7:17151] [11] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(+
> 0x449bf22)[0x7f3800de3f22]
> [c12-7:17151] [12] /share/apps/cig/dealii/dealii-
> 9.0.0/install/deal.II-v9.0.0/lib/libdeal_II.g.so.9.0.0(_
> ZN6dealii11VectorTools11interpolateILi2ELi2ENS_16TrilinosWrapper
> s3MPI11BlockVectorENS_10DoFHandlerEEEvRKNS_7MappingIXT_EXT0_EEERKT2_IXT_
> EXT0_EERKNS_8FunctionIXT0_ENT1_10value_typeEEERSF_RKNS_
> 13ComponentMaskE+0x67)[0x7f3800de47b9]
> [c12-7:17151] [13] /home/billen/AspectProjects/aspect/build/aspect(_
> ZN6aspect9SimulatorILi2EE48set_initial_temperature_and_
> compositional_fieldsEv+0x544)[0x55b228457f18]
> [c12-7:17151] [14] /home/billen/AspectProjects/aspect/build/aspect(_
> ZN6aspect9SimulatorILi2EE3runEv+0xf3d)[0x55b2284028f7]
> [c12-7:17151] [15] /home/billen/AspectProjects/
> aspect/build/aspect(_Z13run_simulatorILi2EEvRKNSt7__
> cxx1112basic_stringIcSt11char_traitsIcESaIcEEEbb+0x117)[0x55b227fda709
> ]
> [c12-7:17151] [16] /home/billen/AspectProjects/aspect/build/aspect(main+
> 0x4df)[0x55b227fb9b21]
> [c12-7:17151] [17] /lib/x86_64-linux-gnu/libc.so.
> 6(__libc_start_main+0xe7)[0x7f37f9f6bb97]
> [c12-7:17151] [18] /home/billen/AspectProjects/aspect/build/aspect(_start+
> 0x2a)[0x55b227e256ea]
> [c12-7:17151] *** End of error message ***
> srun: error: c12-7: task 0: Floating point exception
> srun: error: c12-7: task 1: Exited with exit code 1
> slurmstepd: error: c12-7 [0] pmixp_client_v2.c:203 [_errhandler] mpi/pmix:
> ERROR: Error handler invoked: status = -25: Interrupted system call (4)
> srun: Job step aborted: Waiting up to 32 seconds for job step to finish.
> slurmstepd: error: *** STEP 160.0 ON c12-7 CANCELLED AT
> 2018-06-29T01:20:19 ***
> srun: error: c12-7: tasks 2-7: Killed
> >
>
>
>
> ____________________________________________________________
> Professor of Geophysics
> Earth & Planetary Sciences Dept., UC Davis
> Davis, CA 95616
> 2129 Earth & Physical Sciences Bldg.
> Office Phone: (530) 752-4169
> http://magalibillen.faculty.ucdavis.edu
>
> Currently on Sabbatical at Munich University (LMU)
> Department of Geophysics (PST + 9 hr)
>
> Avoid implicit bias - check before you submit:
> http://www.tomforth.co.uk/genderbias/
> ___________________________________________________________
>
>
> _______________________________________________
> Aspect-devel mailing list
> Aspect-devel at geodynamics.org
> http://lists.geodynamics.org/cgi-bin/mailman/listinfo/aspect-devel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.geodynamics.org/pipermail/aspect-devel/attachments/20180629/79d5ea75/attachment-0001.html>
More information about the Aspect-devel
mailing list