Comments (12)
Would it be possible to include a block of code along the lines of
if (abs(slope).lt.CVMIX_MATH_NEWTON_TOL) then
print*, 'x0 = ', x0
print*, 'root = ', root
print*, 'coeffs(:) = ', coeffs
end if
Just prior to root = root - fun_val/slope
line? Knowing the values of x0
, coeffs(:)
, and root
that are resulting in the divide-by-zero might be helpful in figuring out how to define the bracket if we fall back to a bracketing algorithm.
Relatedly, another ocean parameterization library I'm working on has a logging module that makes it easy to pass status messages / error codes back to the GCM. While it wouldn't fix your problem, I think having something like that (and being able to trace the calls that are resulting in the divide-by-zero, in this case perhaps getting k
values from compute_OBL_depth
) would be extremely useful. So I might implement something like that first, and then try to find a better root-finder afterwards.
I'll try to find some time this weekend or next to work on this, but if you guys make any progress at GFDL please keep me in the loop!
from cvmix-src.
Hi Michael,
Here is some output just prior to the crash. I agree that a return value would be useful in this case. Regarding the bracketing, is there any harm in making this wide (deep)? It could be something that is only done when Newton's method fails.
slope = 0.000000000000000E+000
x0 = -2.05240021255362
root = -2.05240021255362
coeffs(:) = 2.939341553426429E+030 4.296445014121396E+030
2.093375837620301E+030 3.399882447875836E+029
forrtl: error (73): floating divide by zero
Image PC Routine Line Source
MOM6 000000000077D8E5 cvmix_math_mp_cvm 205 cvmix_math.F90
MOM6 0000000000A5E8F4 cvmix_kpp_mp_cvmi 1342 cvmix_kpp.F90
MOM6 000000000073FEEC mom_kpp_mp_kpp_ca 607 MOM_KPP.F90
from cvmix-src.
Michael, would you accept a temporary solution to this problem whereby we do newton's method on the cubic/quadratic and if that fails we fall back to doing it on a linear function? This will involve adding a return value to cvmix_math_cubic_root_find(). It has the advantage of being quick/easy to implement so I can close the MOM6 issue. We can then come up with a more robust solution later (perhaps using brent or similar).
If you're happy with this then I'll send a pull request with the changes.
from cvmix-src.
Falling back to a linear function sounds fine, although adjusting your initial guess might be a good option as well. For this particular case, our initial guess of 0.5*(depth(k)+depth(k+1))
happens to fall directly on a local optimum (or saddle point) of the cubic; I wonder if an initial guess of 0.25*(depth(k)+depth(k+1))
or 0.75*(depth(k)+depth(k+1))
would iterate to the root?
Also, it might be worth investigating why the initial guess is a point where the slope is 0. I'm a little surprised (but I suppose it could be a coincidence), so it would be nice to see depth(k-1:k+1)
and Ri_bulk(k-1:k+1)
(from the call to cvmix_math_poly_interp
at cvmix_kpp.F90:1398
).
from cvmix-src.
depth(k-1:k+1): -2.05240021240362 -2.05240021250362
-2.05240021260362
Ri_bulk(k-1:k+1): -2.739995556420065E-005 -7.835367389899912E-007
0.340014162061937
slope: 0.000000000000000E+000
forrtl: error (65): floating invalid
Image PC Routine Line Source
MOM6 00000000004C5E33 cvmix_math_mp_cvm 200 cvmix_math.F90
MOM6 0000000000576AE8 cvmix_kpp_mp_cvmi 1408 cvmix_kpp.F90
MOM6 00000000004B8893 mom_kpp_mp_kpp_ca 607 MOM_KPP.F90
from cvmix-src.
This is a case where the column has some vanished layers in the interior. The block that needs to catch this situation is:
if (k.eq.size(Ri_bulk)) then
OBL_depth = abs(OBL_limit)
elseif (k.eq.0) then
OBL_depth = abs(zt_cntr(1))
else
if (k.eq.1) then
call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, &
depth(k:k+1), Ri_bulk(k:k+1))
else
call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, &
depth(k:k+1), Ri_bulk(k:k+1), depth(k-1), &
Ri_bulk(k-1))
end if
coeffs(1) = coeffs(1)-CVmix_kpp_params_in%ri_crit
OBL_depth = -cvmix_math_cubic_root_find(coeffs, 0.5_cvmix_r8 * &
(depth(k)+depth(k+1)))
.
.
.
end if
I think something like this would avoid calling the solver with constant depths:
if ( abs( (depth(k+1)-depth(k) ) > 0. ) then ! This could use a tolerance rather than 0., to be robust
OBL_depth = -cvmix_math_cubic_root_find(coeffs, 0.5_cvmix_r8 * &
(depth(k)+depth(k+1)))
else
! Ri_bulk crossed the Ri_crit between k and k+1 but the thickness is vanished
OBL_depth = 0.5_cvmix_r8 * (depth(k)+depth(k+1))
endif
from cvmix-src.
Yeah, I was going to say if depth
doesn't change between levels k-1
and k+1
, then the cubic interpolation is going to be nonsense: you are trying to find cubic polynomial P(x)
such that
P(-2.05240021260362) = -2.739995556420065E-005
P(-2.05240021260362) = -7.835367389899912E-007
P(-2.05240021260362) = 0.340014162061937
and it's computing a 4th constraint that is probably
P'(-2.05240021260362) = inf
I'm not familiar with what you guys are calling "vanishing levels", and wrote most of this code assuming that depth(k) > depth(k+1)
for all k
, so it isn't surprising that have 3 levels with the same depth and different bulk Richardson numbers is causing a problem.
from cvmix-src.
Just to be clear, depth(k) > depth(k+1) does hold, e.g. in this case depth(k-1:k+1) is:
-2.05240021240362
-2.05240021250362
-2.05240021260362
However I see that suitability of the interpolation is probably limited in this case!
To avoid this the tolerance mentioned in @adcroft code suggestion would need to be dependent on the minimum layer thickness.
from cvmix-src.
The highlighted digits make all the difference. I guess better than the tolerance variant above would be a normalization of the local coordinate so that coeffs
are order Ri_bulk
. Something like:
if (k.eq.size(Ri_bulk)) then
OBL_depth = abs(OBL_limit)
elseif (k.eq.0) then
OBL_depth = abs(zt_cntr(1))
else
! dstar is the locally normalized coordinate d* = ( d-depth(k-1) ) / ( depth(k-1) - depth(k+1) )
dstar(1) = 0. ! Corresponds to depth(k-1)
dstar(3) = 1. ! Corresponds to depth(k+1)
if (k>1 .and. depth(k-1) - depth(k+1) > 0.) then
dstar(2) = ( depth(k) - depth(k+1) ) / ( depth(k-1) - depth(k+1) ) ! 0. < dstar(2) < 1.
call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, &
dstar(2:3), Ri_bulk(k:k+1), dstar(1), Ri_bulk(k-1))
else
! For k=1 or vanished separations resort to linear interpolation
dstar(2) = 0.0_cvmix_r8
call cvmix_math_poly_interp(coeffs, CVmix_kpp_params_in%interp_type, &
dstar(2:3), Ri_bulk(k:k+1))
endif
coeffs(1) = coeffs(1)-CVmix_kpp_params_in%ri_crit
OBL_depth = depth( max(1,k-1) ) &
- cvmix_math_cubic_root_find(coeffs, &
0.5_cvmix_r8 * ( dstar(2)+dstar(3) ) ) &
* ( depth( max(1,k-1) ) - depth(k+1) )
.
.
.
end if
from cvmix-src.
Wow, I totally missed the single digit change! I'll use these values of depth
and Ri_bulk
to create a new regression test, and then see how the regression test performs with the normalized depth
values.
from cvmix-src.
Okay, I have a pretty basic test in reg_tests/math/
on mnlevy1981/bugfix/better_interp -- my results aren't identical to yours (probably a combination of compiler differences and the fact that I hard-coded depth
and Ri_bulk
based on your text output above, which might be missing a digit or two). But running my test with gfortran I get output that looks like
Interpolating with depth
----
2.9393415534264454E+030 4.2964450141214149E+030 2.0933758376203048E+030 3.3998824478758336E+029
Polynomial is 0 when depth = -2.0524332880733298
Interpolating with normalized depth
----
1.0199375875839882 -4.7597833756802608 6.7997665835970169 -2.7199066334388067
Polynomial is 0 when depth = -2.0524002125037253
We expect the root to be between -2.05240021260362 and -2.05240021240362; the normalized depth provides a root in this range while the non-normalized depth does not (we only have 5 sig figs in the latter case!). Intel and NAG show a similar degradation in the precision for non-normalized depth, while PGI [16.5] for some reason is the opposite
Interpolating with depth
----
2.9393415534264454E+030 4.2964450141214149E+030 2.0933758376203048E+030
3.3998824478758336E+029
Root when depth = -2.052400212553620
Interpolating with normalized depth
----
1.019937587583988 -4.759783375680261 6.799766583597017
-2.719906633438807
Root when depth = -2.052400212303515
Based on these results, though, I definitely support using the normalized-depth option in cvmix_kpp
. I used the line
norm_depth = (depth - depth(1)) / (depth(3) - depth(1))
to normalize the depth and
root = root*(depth(3) - depth(1)) + depth(1)
to transform back to physical space.
from cvmix-src.
Hi @mnlevy1981 thanks for this. I'll make a pull request based on the work above. Alternatively, if your solution is already complete let me know.
from cvmix-src.
Related Issues (20)
- Update 'd' term in bulk Richardson computation HOT 1
- Min value for unresolved shear HOT 1
- Averaging Nsqr_iface in cvmix_kpp_compute_unresolved_shear HOT 1
- Enhance diffusivity not tied to "MatchBoth" HOT 2
- Remove lavg_N_or_Nsqr HOT 1
- computation of zeta for stable buoyancy forcing and wind stress HOT 6
- Allow convective mixing in the boundary layer? HOT 2
- More efficient tidal mixing HOT 1
- Inconsistency in use of max_nlev
- Another divide-by-zero (caught by CESM in debug mode) HOT 1
- Add continuous integration for testing HOT 1
- The bld/cvmix_setup script uses python 2 HOT 4
- bug in Ekman depth limiter HOT 5
- Error in cvmix_shear.f90
- Error in cvmix_kpp.f90
- TravisCI is failing on netCDF build HOT 1
- License HOT 6
- Default local_mixing_frac value HOT 4
- Optional argument CVmix_conv_params_user used in subroutine. HOT 5
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from cvmix-src.