solver crashes when running 5 days on gx3, 1x1

phil-blain commented 4 years ago

Program received signal SIGFPE, Arithmetic exception. it_nl=0 istep=86 ( 14e pas de temps de la 4e journée)

(gdb) bt
#0  0x000000000093532d in ice_dyn_vp::calc_bvec (nx_block=..., ny_block=..., icellu=..., indxui=..., indxuj=..., stpr=..., cw=..., aiu=..., uarear=..., uocn=..., vocn=..., 
    waterx=..., watery=..., uvel=..., vvel=..., bxfix=..., byfix=..., bx=..., by=..., vrel=...)
    at cicecore/cicedynB/dynamics/ice_dyn_vp.F90:2230
#1  0x00000000008e1dca in ice_dyn_vp::anderson_solver (icellt=..., icellu=..., indxti=..., indxtj=..., indxui=..., indxuj=..., aiu=..., ntot=..., waterx=..., watery=..., 
    bxfix=..., byfix=..., umassdti=..., bvec=..., sol=..., diagvec=..., fpresx=..., fpresy=..., zetad=..., cb=..., halo_info_mask=...)
    at cicecore/cicedynB/dynamics/ice_dyn_vp.F90:807
#2  0x00000000008b3d2e in ice_dyn_vp::imp_solver (dt=...) at cicecore/cicedynB/dynamics/ice_dyn_vp.F90:447
#3  0x00000000012e0a82 in ice_step_mod::step_dyn_horiz (dt=...) at cicecore/cicedynB/general/ice_step_mod.F90:627
#4  0x00000000004149aa in cice_runmod::ice_step () at cicecore/drivers/cice/CICE_RunMod.F90:240
#5  0x0000000000413a73 in cice_runmod::cice_run () at cicecore/drivers/cice/CICE_RunMod.F90:77
(gdb) p i
$35 = 65
(gdb) p j
$36 = 95

stPr(i,j), is NaN, which is weird because it's supposed to be correctly initialized in calc_zeta_Pr... this grid point is a point next to the southern coast of Alaska (if those i j are the same as the ones Ncview gives)

also very weird that this did not happen when running the 5 year test on daley (but that was 8x1, so a different decomposition..., so maybe not that weird)

@JFLemieux73

phil-blain commented 4 years ago

See https://github.com/phil-blain/CICE/issues/15#issuecomment-614119749

JFLemieux73 commented 4 years ago

C'est la cause du crash?

phil-blain commented 4 years ago

oui. À ce point dans le calcul, il utilise une valeur du tableau stPr qui n'est pas initialisée, et je compile avec -init=snans,arrays donc il abort quand il tombe sur le NaN... mais je ne comprends pas vraiment pourquoi encore, normalement elle devrait être initilialisée.

phil-blain commented 4 years ago

le tableau n'est pas initialisé à cet endroit parce que l'initialisation de stPr se fait dans ¢alc_zeta_Pr avec un boucle sur icellt, mais le tableau est utilisé dans calc_bvec dans une boucle sur icellu. Mon raisonnement ici https://github.com/phil-blain/CICE/issues/15#issuecomment-614170344 est erroné parce que les deux masques ne sont pas définis de la même façon:

iceumask est défini dans ice_dyn_shared::dyn_prep2 à partir de umask, aiu et umass
icetmask est défini dans ice_dyn_shared::dyn_prep1

donc:

ça veut dire qu'à ce point de la simulation il y a une différence entre les points de glace sur la grille u et sur la grille v.
est-ce que c'est ok que le tableau soit initialisé sur icellt et utilisé sur icellu ? (est-ce que une des deux boucles devrait être changée?)

JFLemieux73 commented 4 years ago

Tu veux dire grille u et grille t?

Faudrait que je checke. Il y a parfois des trucs de la grille t (les quatre autour) qui sont utilisés pour les calculs au point u. C'est peut-être ça qui se passe.

jf

JFLemieux73 commented 4 years ago

Selon moi ca crash dans le marginal ice zone. Il y a très peu de glace (Aice tend vers zero). Il y a un calcul au point u (Aice > puny) qui a besoin des 4 voisins au point t. Un des point t est un nan car Aice < puny. La même chose peut se passer avec le evp mais ca doit être initialisé avec des zéros au lieu des nan. Je pense que ça marcherait en l'initialisant avec des zéros. D'un point de vue physique c'est ok que les terms de rhéologie soient à zéro vue le peu de glace.

phil-blain commented 4 years ago

Merci JF, je pense que tu as raison. J'ai regardé le code du EVP, et effectivement le tableau str utilisé dans ice_dyn_shared::stepu, qui intègre l'équation de momentum, https://github.com/phil-blain/CICE/blob/8bf95d102c6779f227e431026768e7827e79abdf/cicecore/cicedynB/dynamics/ice_dyn_shared.F90#L722-L726 est initialisé à zéro dans ice_dyn_evp::stress: https://github.com/phil-blain/CICE/blob/8bf95d102c6779f227e431026768e7827e79abdf/cicecore/cicedynB/dynamics/ice_dyn_evp.F90#L626-L629

(avant de calculer explicitement les stress pour les cellules où il y a de la glace).

Donc je vais faire la même chose pour notre solveur, initialiser stPr à zéro. Merci d'avoir confirmé que ça fait du sens niveau physique.

phil-blain commented 4 years ago

@JFLemieux73 j'imagine que ça fait aussi du sens que zetaD soit initialisé à zéro ?

JFLemieux73 commented 4 years ago

oui car zeta=ice_strength / (2 Delta*) ...dans ces zones avec très peu de glace, ice_strength tend vers zéro.

phil-blain commented 4 years ago

Avec l'initialisation de zetaD et stPr, la code complète avec succès. Je reprends donc les tests MPI.

phil-blain / CICE

solver crashes when running 5 days on gx3, 1x1 #34