From 9cfd382c6a340c4ad56db5c90f82aec6158a31c5 Mon Sep 17 00:00:00 2001 From: APiccolo89 <85160533+APiccolo89@users.noreply.github.com> Date: Thu, 8 Feb 2024 14:06:51 +0100 Subject: [PATCH 1/3] Bug mode restart 1.] This bug creates problem only under these circumstances: if someone activates the diffusive temperature at the beginning of the simulation/ if someone is activating the lithostatic pressure flag. -> If you the simulation is restarted, the temperature diffusion step is repeat at the beginning of the restarting simulation, as well as for the lithostatic pressure. => Since the lithostatic pressure is substituting the actual pressure field for the computation, this invalidates a bit the restarted simulation. Same reasoning can be applied with the temperature diffusive timestep. --- src/JacRes.cpp | 10 ++++++---- src/JacRes.h | 4 ++-- src/LaMEMLib.cpp | 8 +++++--- 3 files changed, 13 insertions(+), 9 deletions(-) diff --git a/src/JacRes.cpp b/src/JacRes.cpp index c10d6202..7357125e 100644 --- a/src/JacRes.cpp +++ b/src/JacRes.cpp @@ -1930,9 +1930,11 @@ PetscErrorCode JacResCopyPres(JacRes *jr, Vec x) PetscFunctionReturn(0); } //--------------------------------------------------------------------------- -PetscErrorCode JacResInitPres(JacRes *jr) +PetscErrorCode JacResInitPres(JacRes *jr,TSSol *ts) + { FDSTAG *fs; + BCCtx *bc; SolVarCell *svCell; const PetscScalar *p; @@ -1949,7 +1951,7 @@ PetscErrorCode JacResInitPres(JacRes *jr) fixPhase = bc->fixPhase; // check activation - if(!bc->initPres) PetscFunctionReturn(0); + if(!bc->initPres || ts->istep>0) PetscFunctionReturn(0); // get grid coordinate bounds in z-direction ierr = FDSTAGGetGlobalBox(fs, NULL, NULL, &bz, NULL, NULL, &ez); CHKERRQ(ierr); @@ -1999,7 +2001,7 @@ PetscErrorCode JacResInitPres(JacRes *jr) PetscFunctionReturn(0); } //--------------------------------------------------------------------------- -PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx) +PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx,TSSol *ts) { FDSTAG *fs; SolVarCell *svCell; @@ -2015,7 +2017,7 @@ PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx) PetscFunctionBeginUser; // check activation - if(!jr->ctrl.initLithPres) PetscFunctionReturn(0); + if(!jr->ctrl.initLithPres || ts->istep >0) PetscFunctionReturn(0); // print PrintStart(&t, "Initializing pressure with lithostatic pressure", NULL); diff --git a/src/JacRes.h b/src/JacRes.h index 994b3f76..833c0113 100644 --- a/src/JacRes.h +++ b/src/JacRes.h @@ -313,10 +313,10 @@ PetscErrorCode JacResCopyVel(JacRes *jr, Vec x); PetscErrorCode JacResCopyPres(JacRes *jr, Vec x); // initialize pressure -PetscErrorCode JacResInitPres(JacRes *jr); +PetscErrorCode JacResInitPres(JacRes *jr,TSSol *ts); // initialize pressure to lithostatic pressure -PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx); +PetscErrorCode JacResInitLithPres(JacRes *jr, AdvCtx *actx, TSSol *ts); // copy residuals from local to global vectors, enforce boundary constraints PetscErrorCode JacResCopyRes(JacRes *jr, Vec f); diff --git a/src/LaMEMLib.cpp b/src/LaMEMLib.cpp index d30e6636..e4fd7aee 100644 --- a/src/LaMEMLib.cpp +++ b/src/LaMEMLib.cpp @@ -822,10 +822,10 @@ PetscErrorCode LaMEMLibInitGuess(LaMEMLib *lm, SNES snes) ierr = LaMEMLibDiffuseTemp(lm); CHKERRQ(ierr); // initialize pressure - ierr = JacResInitPres(&lm->jr); CHKERRQ(ierr); + ierr = JacResInitPres(&lm->jr,&lm->ts); CHKERRQ(ierr); // lithostatic pressure initializtaion - ierr = JacResInitLithPres(&lm->jr, &lm->actx); CHKERRQ(ierr); + ierr = JacResInitLithPres(&lm->jr, &lm->actx, &lm->ts); CHKERRQ(ierr); // compute inverse elastic parameters (dependent on dt) ierr = JacResGetI2Gdt(&lm->jr); CHKERRQ(ierr); @@ -864,6 +864,7 @@ PetscErrorCode LaMEMLibInitGuess(LaMEMLib *lm, SNES snes) PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm) { JacRes *jr; + TSSol *ts; Controls *ctrl; AdvCtx *actx; PetscLogDouble t; @@ -874,12 +875,13 @@ PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm) PetscFunctionBeginUser; // access context + ts = &lm->ts; jr = &lm->jr; ctrl = &jr->ctrl; actx = &lm->actx; // check for infinite diffusion - if (ctrl->actTemp && ctrl->actSteadyTemp) + if (ctrl->actTemp && ctrl->actSteadyTemp && ts->istep<=1) { PrintStart(&t,"Computing steady-state temperature distribution", NULL); From ee0df5aa3acd307538c42d69ff8935adb5b677ba Mon Sep 17 00:00:00 2001 From: APiccolo89 <85160533+APiccolo89@users.noreply.github.com> Date: Thu, 8 Feb 2024 14:25:11 +0100 Subject: [PATCH 2/3] Better istep ==0 Just checking now if it works --- src/LaMEMLib.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/LaMEMLib.cpp b/src/LaMEMLib.cpp index e4fd7aee..8b87eca7 100644 --- a/src/LaMEMLib.cpp +++ b/src/LaMEMLib.cpp @@ -629,7 +629,7 @@ PetscErrorCode LaMEMLibSolve(LaMEMLib *lm, void *param, PetscLogStage stages[4]) PetscCall(PetscLogStagePush(stages[0])); /* Start profiling stage*/ ierr = LaMEMLibInitGuess(lm, snes); CHKERRQ(ierr); - + PetscCall(PetscLogStagePop()); /* Stop profiling stage*/ if (param) @@ -881,7 +881,7 @@ PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm) actx = &lm->actx; // check for infinite diffusion - if (ctrl->actTemp && ctrl->actSteadyTemp && ts->istep<=1) + if (ctrl->actTemp && ctrl->actSteadyTemp && ts->istep==0) { PrintStart(&t,"Computing steady-state temperature distribution", NULL); From 3042c5448118b48c5b426f4f1e7fdf2b130a7c5d Mon Sep 17 00:00:00 2001 From: APiccolo89 <85160533+APiccolo89@users.noreply.github.com> Date: Thu, 8 Feb 2024 14:46:13 +0100 Subject: [PATCH 3/3] Solved. Solved. --- src/LaMEMLib.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/LaMEMLib.cpp b/src/LaMEMLib.cpp index 8b87eca7..589dd706 100644 --- a/src/LaMEMLib.cpp +++ b/src/LaMEMLib.cpp @@ -905,7 +905,7 @@ PetscErrorCode LaMEMLibDiffuseTemp(LaMEMLib *lm) } // check for additional limited diffusion - if (ctrl->actTemp && ctrl->steadyTempStep) + if (ctrl->actTemp && ctrl->steadyTempStep && ts->istep==0) { PrintStart(&t,"Diffusing temperature", NULL);