v1.1.0 release

.gitignore

@@ -2,8 +2,11 @@ debug
 build
 tests/*/*.F90
 .VSCodeCounter
-.vscode
+.vscode*
 ford_doc.md
 *mod
 *.json
 /doc/
+.ssh
+.gitconfig
+*bash_history

CHANGELOG.md

@@ -1,15 +1,29 @@
 # CHANGELOG
 
-## v1.1.0, 2023-06-28
+## v1.1.0, 2024-02-02
 
 - Solver and integrator improvements
+- Debug mode bug fixes
+- Multilinear interpolation support
+- Support for some CRM features with non-default electron species ID
+- Default support for release build assertions during startup
 
 ### New Features
 
 - The BDE integrator internal controller now attempts to reduce the number of substeps back to 1 after every 50 successful integrations
-- Command line PETSc support for setting up the KSP solver object has been implemented. The command line ksp_type takes precedent over the integrator 
-options. This now allows for command line customization of the KSP object. Note that tolerances are still set from config options. 
+- The BDE integrator now also has a hard maximum of attempted restarts regardless of whether consolidation happens. this avoids crashes due to allocating too much memory for the timestep size buffer
+- The BDE integrator now prints the approximate convergence bottleneck variable (feature not completely reliable)
+- Command line PETSc support for setting up the KSP solver object has been implemented. The command line ksp_type takes precedent over the integrator options. This now allows for command line customization of the KSP object. Note that tolerances are still set from config options. 
+- New derivation type nDLinInterpDerivation - linear interpolation on n-dimensional data
+- ModelboundCRMData can now have an associated non-default electron species ID (useful for multiple flux tube models). Note that Janev transitions still have hardcoded values. 
+- Added the option to reset the value of the time variable upon loading from restart files
+- Assertions are now run in release mode during simulation startup
 
+### Bug Fixes
+
+- Fixed a number of bugs in the code where divide-by-zero FPEs would be raised in debug mode. Some still remain.
+- Fixed a bug with gfortran-11.4 which caused bound changes in an MPI buffer
+- Some tests were failing when run on an M2 Mac due to FPE differences. This has been fixed. 
 
 ## v1.0.0, 2023-06-21
 

CMakeLists.txt

@@ -17,7 +17,7 @@
 
 CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
 
-PROJECT(ReMKiT1D VERSION 1.0.0
+PROJECT(ReMKiT1D VERSION 1.1.0
   DESCRIPTION "A Fortran framework for 1D multifluid and kinetic simulations geared towards Scrape-Off Layer plasmas"
   LANGUAGES Fortran)
 

codecov.yml

@@ -0,0 +1,16 @@
+coverage:
+  round: up
+  range: 50..80
+
+  status:
+    project: # settings affecting project coverage
+      default:
+        enabled: yes
+        target: auto # auto % coverage target
+        threshold: 5%  # allow for 5% reduction of coverage without failing
+
+    # do not run coverage on patch nor changes
+    patch:
+      default:
+        enabled: false
+    changes: no

src/modules/basic_integrators/explicit_rk_integrator.f90

@@ -18,7 +18,7 @@ module explicit_rk_integrator_class
 
     use data_kinds                         ,only: rk ,ik
     use god_objects                        ,only: Object
-    use runtime_constants                  ,only: debugging, assertions
+    use runtime_constants                  ,only: debugging, assertions, assertionLvl
     use assertion_utility                  ,only: assert, assertIdentical, assertPure
     use modeller_surrogate_class           ,only: ModellerSurrogate
     use variable_container_class           ,only: VariableContainer

src/modules/basic_integrators/explicit_rk_integrator_procedures.f90

@@ -39,7 +39,7 @@ module subroutine initRKIntegrator(this,modelList,termGroups,order,tableau,evolv
 
     integer(ik)                                         :: i,j
 
-    if (assertions) then
+    if (assertions .or. assertionLvl >= 0) then
 
         if (present(termGroups)) call assert(present(modelList),"Term groups object passed to RK integrator constructor without&
                                                                & model list")
@@ -219,7 +219,7 @@ module subroutine integrateRK(this,manipulatedModeller,outputVars,inputVars)
                 if (.not. mixedGroup(i)%entry(j)) then
                     varIndices(i)%entry(j) = inputVars%getVarIndex(manipulatedModeller&
                                                               %getEvolvedVarInTermGroup(termGroups(i)%entry(j),modelIndices(i)))
-                    if (assertions) call assert(.not.inputVars%isStationary(manipulatedModeller&
+                    call assert(.not.inputVars%isStationary(manipulatedModeller&
                     %getEvolvedVarInTermGroup(termGroups(i)%entry(j),modelIndices(i))),&
                     "Explicit RK integrator detected stationary variable among evolved variables - this is unsupported")
                 end if

src/modules/basic_integrators/implicit_PicardBDE_integrator.f90

@@ -18,7 +18,7 @@ module implicit_PicardBDE_integrator_class
 
     use data_kinds                         ,only: rk ,ik
     use god_objects                        ,only: Object
-    use runtime_constants                  ,only: debugging, assertions
+    use runtime_constants                  ,only: debugging, assertions, assertionLvl
     use assertion_utility                  ,only: assert, assertIdentical, assertPure
     use modeller_surrogate_class           ,only: ModellerSurrogate
     use variable_container_class           ,only: VariableContainer
@@ -42,6 +42,9 @@ module implicit_PicardBDE_integrator_class
         integer(ik) :: minNonlinIters = 5 !! Number of nonlinear iterations below which the number of substeps gets reduced
 
         integer(ik) :: maxRestarts = 3 !! Maximum number of consecutive solver restart attempts before critical failure is announced
+
+        integer(ik) :: hardMaxRestarts = 10 !! Maximum number of consecutive solver restart attempts before critical failure is announced regardless of whether consolidation happened or not
+
         integer(ik) :: restartCount = 0 !! Counter for consecutive number of solver restars 
 
         integer(ik) :: consolidationInterval = 50 !! How many integration calls before currentNumSubsteps is again reduced to 1

src/modules/basic_integrators/implicit_PicardBDE_integrator_procedures.f90

@@ -45,7 +45,7 @@ module subroutine initBDEIntegrator(this,indexingObj,procRank,nonlinTol,absTol,m
 
     integer(ik) ,allocatable ,dimension(:) :: procDoFs
 
-    if (assertions) then
+    if (assertions .or. assertionLvl >= 0) then
         call assert(indexingObj%isDefined(),"Undefined indexing object passed to BDE integrator constructor")
         if (present(termGroups)) call assert(present(modelList),"Term groups object passed to BDE integrator constructor without&
                                                                & model list")
@@ -146,7 +146,10 @@ module subroutine integrateBDE(this,manipulatedModeller,outputVars,inputVars)
             allocate(dt(1))
             dt = fullTimestep
         end if
-        if (this%internalControlOpts%restartCount > this%internalControlOpts%maxRestarts) error stop "Max BDE restarts reached"
+        if (this%internalControlOpts%restartCount > this%internalControlOpts%hardMaxRestarts) &
+            error stop "Max BDE restarts reached - hard maximum"
+        if (this%internalControlOpts%restartCount > this%internalControlOpts%maxRestarts .and. &
+            this%internalControlOpts%stepsSinceLastConsolidation > 1) error stop "Max BDE restarts reached"
         call tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,dt,solveSuccess)
         if (solveSuccess) then
             this%internalControlOpts%restartCount = 0
@@ -163,7 +166,7 @@ module subroutine integrateBDE(this,manipulatedModeller,outputVars,inputVars)
 
 end subroutine integrateBDE
 !-----------------------------------------------------------------------------------------------------------------------------------
-function checkConvergence(oldVars,newVars,indicesToCheck,nonlinTol,absTol,use2Norm) result(conv)
+function checkConvergence(oldVars,newVars,indicesToCheck,nonlinTol,absTol,use2Norm,convergenceCounter) result(conv)
     !! Checks whether all variables determined by indicesToCheck have converged based on a given nonlinear tolerance
 
     type(RealArray) ,dimension(:) ,intent(in) :: oldVars !! Previous variable values
@@ -172,7 +175,7 @@ function checkConvergence(oldVars,newVars,indicesToCheck,nonlinTol,absTol,use2No
     real(rk)                      ,intent(in) :: nonlinTol !! Relative convergence tolerances for each variable 
     real(rk)                      ,intent(in) :: absTol !! Absolute tolerance in epsilon units for each variable
     logical                       ,intent(in) :: use2Norm !! True if this should use 2-norm instead of local inf-norm
-
+    integer(ik)     ,dimension(:) ,intent(inout) :: convergenceCounter !! Incremented if the variable with a given index has not yet converged
 
     logical :: conv 
     logical ,allocatable ,dimension(:) :: varConverged
@@ -188,15 +191,29 @@ function checkConvergence(oldVars,newVars,indicesToCheck,nonlinTol,absTol,use2No
         nonHaloLen = size(newVars(ind)%entry) - haloDataChunkSize
 
         if (use2Norm) then
-            relError = norm2(oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen))&
-                    /norm2(oldVars(ind)%entry(1:nonHaloLen))
+
             absError = norm2(oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen))
+
+            if (norm2(oldVars(ind)%entry(1:nonHaloLen))>epsilon(absError)*absTol) then
+                relError = norm2(oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen))&
+                        /norm2(oldVars(ind)%entry(1:nonHaloLen))
+            else
+                relError = absError
+            end if
         else
-            relError = maxval(abs((oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen)))&
-                    /abs(oldVars(ind)%entry(1:nonHaloLen)))
+
             absError = maxval(abs((oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen))))
+
+            if (all(abs(oldVars(ind)%entry(1:nonHaloLen)) > epsilon(absError)*absTol)) then 
+                relError = maxval(abs((oldVars(ind)%entry(1:nonHaloLen)-newVars(ind)%entry(1:nonHaloLen)))&
+                    /abs(oldVars(ind)%entry(1:nonHaloLen)))
+                else
+                    relError = absError
+            end if
+
         end if
         varConverged(i) = relError < nonlinTol .or. absError < epsilon(absError)*absTol
+        if (.not. varConverged(i)) convergenceCounter(i) = convergenceCounter(i) + 1
     end do
 
     conv = all(varConverged)
@@ -300,6 +317,8 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
     real(rk)        ,allocatable ,dimension(:) :: implicitVectorInit
 
     real (rk) :: startingTime 
+
+    integer(ik)    ,allocatable ,dimension(:)  :: convergenceCounter 
 
     termGroups = this%getTermGroups()
     modelIndices = this%getModelIndices()
@@ -309,6 +328,8 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
     if (nonTrivialModelDataUpdate) modelDataUpdateRules = this%getModelDataUpdateRules()
     nonTrivialConvergenceCheck = allocated(this%convergenceTestVars) 
 
+    if (nonTrivialConvergenceCheck) allocate(convergenceCounter(size(this%convergenceTestVars)))
+
     !Check if communication needed
     commNeeded = this%isCommunicationNeeded()
     if (commNeeded) commData = this%getCommunicationData()
@@ -330,7 +351,7 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
             if (inputVars%isVarNameRegistered("time")) &
             this%buffer%variables(timeVarIndex)%entry(1) = this%buffer%variables(timeVarIndex)%entry(1) + dt(i)
             tolReached = .false.
-
+            if (nonTrivialConvergenceCheck) convergenceCounter = 1
 
             call this%buffer%copyImplicitVarsToVec(implicitVectorInit,ignoreStationary=.true.)
             do nonlinIter = 1, this%maxIterations
@@ -371,7 +392,7 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
                         this%internalControlOpts%currentNumSubsteps = this%internalControlOpts%currentNumSubsteps &
                                                                     * this%internalControlOpts%stepMultiplier
                         solveSuccess = .false.
-                        this%internalControlOpts%restartCount = this%internalControlOpts%restartCount
+                        this%internalControlOpts%restartCount = this%internalControlOpts%restartCount + 1
                         return
                     else
                         call printNamedValue(this%integratorName//": PETScConvergedReason:",convReason)
@@ -391,7 +412,7 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
                 if (nonTrivialConvergenceCheck) then
                     locConverged = &
                     checkConvergence(oldBufferVals,this%buffer%variables,this%convergenceTestVars,&
-                                    this%nonlinTol,this%absTol,this%use2Norm)
+                                    this%nonlinTol,this%absTol,this%use2Norm,convergenceCounter)
                 else
                     locConverged = (norm2(this%implicitVectorOld-this%implicitVectorNew)/norm2(this%implicitVectorOld)) &
                                 < this%nonlinTol
@@ -413,6 +434,14 @@ subroutine tryIntegrate(this,manipulatedModeller,outputVars,inputVars,numSteps,d
                 end if
             end if
 
+            if (nonTrivialConvergenceCheck) then
+                do j = 1, size(convergenceCounter)
+                    if (convergenceCounter(j) == nonlinIter) &
+                        call printMessage(this%integratorName//": convergence bottleneck: "&
+                                    // inputVars%getVarName(this%convergenceTestVars(j)),.true.)
+                end do
+            end if
+
             if (solveSuccess) then
                 this%totNumIters = this%totNumIters + nonlinIter
                 this%timesCalled = this%timesCalled + 1

src/modules/basic_integrators/simple_timestep_controller.f90

@@ -17,7 +17,7 @@ module simple_timestep_controller_class
     !! Houses simple timestep controller which uses the min (max) value of a product of variables, scaled to corresponding powers
 
     use data_kinds                         ,only: rk ,ik
-    use runtime_constants                  ,only: debugging, assertions
+    use runtime_constants                  ,only: debugging, assertions, assertionLvl
     use god_objects                        ,only: Object
     use assertion_utility                  ,only: assert, assertIdentical, assertPure
     use variable_container_class           ,only: VariableContainer

src/modules/basic_integrators/simple_timestep_controller_procedures.f90

@@ -36,7 +36,7 @@ module subroutine initSimpleTimestepController(this,mpiCont,varCont,reqVars,reqV
 
     integer(ik) :: i
 
-    if (assertions) then 
+    if (assertions .or. assertionLvl >= 0) then 
         call assert(mpiCont%isDefined(),"Undefined MPI controller passed to initSimppleTimestepController")
         call assert(varCont%isDefined(),"Undefined reference variable container passed to initSimppleTimestepController")
         call assert(size(reqVars) == size(reqVarPowers), "reqVars and reqVarPowers passed to initSimpleTimestepController &