Add convergence to the monitor

src/.depends

@@ -197,7 +197,7 @@ fluid/stress_formulation/pnpn_res_stress_fctry.lo : fluid/stress_formulation/pnp
 fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.lo : fluid/stress_formulation/bcknd/cpu/pnpn_res_stress_cpu.f90 math/math.lo sem/space.lo config/num_types.lo mesh/mesh.lo field/scratch_registry.lo fluid/pnpn_res.lo bc/facet_normal.lo sem/coef.lo math/ax.lo field/field.lo math/operators.lo gs/gather_scatter.lo 
 fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.lo : fluid/stress_formulation/bcknd/device/pnpn_res_stress_device.F90 math/bcknd/device/device_math.lo math/bcknd/device/device_mathops.lo sem/space.lo config/num_types.lo mesh/mesh.lo field/scratch_registry.lo fluid/pnpn_res.lo bc/facet_normal.lo sem/coef.lo math/ax.lo field/field.lo math/operators.lo common/utils.lo gs/gather_scatter.lo 
 fluid/fluid_scheme.lo : fluid/fluid_scheme.f90 common/gradient_jump_penalty.lo bc/shear_stress.lo bc/wall_model_bc.lo math/field_math.lo common/time_step_controller.lo field/field_series.lo common/utils.lo common/user_intf.lo field/scratch_registry.lo common/json_utils.lo field/field_registry.lo common/log.lo math/operators.lo time_schemes/time_scheme_controller.lo math/bcknd/device/device_math.lo math/math.lo mesh/mesh.lo bc/bc.lo fluid/fluid_stats.lo krylov/precon.lo krylov/pc_hsmg.lo krylov/bcknd/device/pc_jacobi_device.lo krylov/bcknd/sx/pc_jacobi_sx.lo krylov/bcknd/cpu/pc_jacobi.lo bc/field_dirichlet_vector.lo bc/field_dirichlet.lo bc/symmetry.lo bc/dong_outflow.lo bc/dirichlet.lo bc/blasius.lo bc/usr_inflow.lo bc/inflow.lo bc/wall.lo sem/coef.lo krylov/krylov.lo sem/dofmap.lo sem/space.lo field/field.lo fluid/fluid_source_term.lo comm/comm.lo config/num_types.lo fluid/mean_flow.lo common/checkpoint.lo config/neko_config.lo fluid/mean_sqr_flow.lo gs/gather_scatter.lo 
-fluid/fluid_aux.lo : fluid/fluid_aux.f90 krylov/krylov.lo config/num_types.lo common/log.lo 
+fluid/fluid_aux.lo : fluid/fluid_aux.f90 common/utils.lo krylov/krylov.lo config/num_types.lo common/log.lo 
 fluid/fluid_pnpn.lo : fluid/fluid_pnpn.f90 math/operators.lo math/field_math.lo common/utils.lo bc/bc.lo math/mathops.lo math/math.lo config/neko_config.lo gs/gs_ops.lo common/time_step_controller.lo common/user_intf.lo mesh/mesh.lo comm/comm.lo bc/non_normal.lo bc/facet_normal.lo bc/dirichlet.lo field/field.lo math/ax.lo common/json_utils.lo common/profiler.lo fluid/advection.lo device/device.lo common/projection.lo time_schemes/time_scheme_controller.lo fluid/fluid_aux.lo math/bcknd/device/device_mathops.lo fluid/fluid_scheme.lo fluid/fluid_volflow.lo common/rhs_maker.lo fluid/pnpn_res.lo krylov/krylov.lo config/num_types.lo 
 fluid/fluid_fctry.lo : fluid/fluid_fctry.f90 common/utils.lo fluid/fluid_pnpn.lo fluid/fluid_scheme.lo 
 fluid/fluid_volflow.lo : fluid/fluid_volflow.f90 math/ax.lo bc/bc.lo field/scratch_registry.lo common/json_utils.lo gs/gather_scatter.lo math/bcknd/device/device_mathops.lo math/bcknd/device/device_math.lo config/neko_config.lo comm/comm.lo math/math.lo time_schemes/time_scheme_controller.lo sem/coef.lo field/field.lo sem/dofmap.lo krylov/precon.lo krylov/krylov.lo math/mathops.lo config/num_types.lo math/operators.lo 

src/fluid/fluid_aux.f90

@@ -4,6 +4,7 @@ module fluid_aux
   use num_types, only : rp
   use krylov, only : ksp_monitor_t
   use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
+  use utils, only : neko_error, neko_warning
   implicit none
   private
 
@@ -13,14 +14,14 @@ module fluid_aux
 
   !> Prints for prs, velx, vely, velz the following:
   !! Number of iterations, start residual, end residual
-  subroutine fluid_step_info(step, t, dt, ksp_results)
+  subroutine fluid_step_info(step, t, dt, ksp_results, strict_convergence)
     type(ksp_monitor_t), intent(in) :: ksp_results(4)
     integer, intent(in) :: step
     real(kind=rp), intent(in) :: t, dt
+    logical, intent(in) :: strict_convergence
     character(len=LOG_SIZE) :: log_buf
     integer :: i
 
-
     call neko_log%message('Pressure')
 
     write(log_buf, '(A,A,A)') 'Iterations:   ',&
@@ -54,11 +55,30 @@ subroutine fluid_step_info(step, t, dt, ksp_results)
          ksp_results(4)%res_start, ksp_results(4)%res_final
     call neko_log%message(log_buf)
 
-    ! Check for divergence
+    ! Check for convergence
     do i = 1, 4
        if (ieee_is_nan(ksp_results(i)%res_final)) then
-          call neko_log%error("Fluid solver diverged")
-          stop
+          call neko_error("Fluid solver diverged")
+       end if
+
+       if (.not. ksp_results(i)%converged) then
+          log_buf = 'Fluid solver did not converge for '
+          select case(i)
+            case(1)
+             log_buf = trim(log_buf) // 'pressure'
+            case(2)
+             log_buf = trim(log_buf) // 'x-velocity'
+            case(3)
+             log_buf = trim(log_buf) // 'y-velocity'
+            case(4)
+             log_buf = trim(log_buf) // 'z-velocity'
+          end select
+
+          if (strict_convergence) then
+             call neko_error(log_buf)
+          else
+             call neko_warning(log_buf)
+          end if
        end if
     end do
 

src/fluid/fluid_pnpn.f90

@@ -822,7 +822,7 @@ subroutine fluid_pnpn_step(this, t, tstep, dt, ext_bdf, dt_controller)
               this%ksp_vel%max_iter)
       end if
 
-      call fluid_step_info(tstep, t, dt, ksp_results)
+      call fluid_step_info(tstep, t, dt, ksp_results, this%strict_convergence)
 
       call this%scratch%relinquish_field(temp_indices)
 

src/fluid/fluid_scheme.f90

@@ -104,14 +104,18 @@ module fluid_scheme
      type(field_t), pointer :: f_y => null()
      !> Z-component of the right-hand side.
      type(field_t), pointer :: f_z => null()
-     class(ksp_t), allocatable :: ksp_vel     !< Krylov solver for velocity
-     class(ksp_t), allocatable :: ksp_prs     !< Krylov solver for pressure
-     class(pc_t), allocatable :: pc_vel        !< Velocity Preconditioner
-     class(pc_t), allocatable :: pc_prs        !< Pressure Preconditioner
+
+     ! Krylov solvers and settings
+     class(ksp_t), allocatable :: ksp_vel  !< Krylov solver for velocity
+     class(ksp_t), allocatable :: ksp_prs  !< Krylov solver for pressure
+     class(pc_t), allocatable :: pc_vel    !< Velocity Preconditioner
+     class(pc_t), allocatable :: pc_prs    !< Velocity Preconditioner
      integer :: vel_projection_dim         !< Size of the projection space for ksp_vel
      integer :: pr_projection_dim          !< Size of the projection space for ksp_pr
      integer :: vel_projection_activ_step  !< Steps to activate projection for ksp_vel
      integer :: pr_projection_activ_step   !< Steps to activate projection for ksp_pr
+     logical :: strict_convergence         !< Strict convergence for the velocity solver
+
      type(no_slip_wall_t) :: bc_wall           !< No-slip wall for velocity
      class(bc_t), allocatable :: bc_inflow !< Dirichlet inflow for velocity
      type(wall_model_bc_t) :: bc_wallmodel !< Wall model boundary condition
@@ -675,6 +679,10 @@ subroutine fluid_scheme_init_common(this, msh, lx, params, scheme, user, &
        call neko_log%end_section()
     end if
 
+    ! Strict convergence for the velocity solver
+    call json_get_or_default(params, 'case.fluid.strict_convergence', &
+         this%strict_convergence, .false.)
+
     ! Assign velocity fields
     call neko_field_registry%add_field(this%dm_Xh, 'u')
     call neko_field_registry%add_field(this%dm_Xh, 'v')

src/krylov/bcknd/cpu/bicgstab.f90

@@ -244,6 +244,7 @@ function bicgstab_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
   end function bicgstab_solve
 
   !> Standard BiCGSTAB coupled solve

src/krylov/bcknd/cpu/cacg.f90

@@ -328,6 +328,7 @@ function cacg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resul
       call this%monitor_stop()
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
 
     end associate
 

src/krylov/bcknd/cpu/cg.f90

@@ -233,7 +233,7 @@ function cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_results
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
-
+    ksp_results%converged = this%is_converged(iter, rnorm)
   end function cg_solve
 
   subroutine second_cg_part(rtr, r, mult, w, alpha, n)

src/krylov/bcknd/cpu/cg_coupled.f90

@@ -325,6 +325,7 @@ function cg_cpld_solve(this, Ax, x, y, z, fx, fy, fz, &
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
   end function cg_cpld_solve
 
 end module cg_cpld

src/krylov/bcknd/cpu/cheby.f90

@@ -240,6 +240,7 @@ function cheby_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
       rnorm = sqrt(rtr) * norm_fac
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
     end associate
   end function cheby_solve
 

src/krylov/bcknd/cpu/gmres.f90

@@ -364,6 +364,7 @@ function gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
 
   end function gmres_solve
 

src/krylov/bcknd/cpu/pipecg.f90

@@ -362,6 +362,7 @@ function pipecg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_res
       call this%monitor_stop()
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
 
     end associate
 

src/krylov/bcknd/device/cg_device.f90

@@ -230,6 +230,7 @@ function cg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
 
   end function cg_device_solve
 

src/krylov/bcknd/device/cheby_device.F90

@@ -285,8 +285,11 @@ function cheby_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
       call device_sub2(r_d, w_d, n)
       rtr = device_glsc3(r_d, coef%mult_d, r_d, n)
       rnorm = sqrt(rtr) * norm_fac
+
+
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
     end associate
   end function cheby_device_solve
 

src/krylov/bcknd/device/fusedcg_cpld_device.F90

@@ -635,6 +635,7 @@ function fusedcg_cpld_device_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
       call this%monitor_stop()
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
 
     end associate
 
@@ -659,6 +660,7 @@ function fusedcg_cpld_device_solve(this, Ax, x, f, n, coef, blst, &
 
     ksp_results%res_final = 0.0
     ksp_results%iter = 0
+    ksp_results%converged = .false.
 
   end function fusedcg_cpld_device_solve
 

src/krylov/bcknd/device/fusedcg_device.F90

@@ -394,6 +394,7 @@ function fusedcg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result
       call this%monitor_stop()
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
 
     end associate
 

src/krylov/bcknd/device/gmres_device.F90

@@ -479,6 +479,7 @@ function gmres_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) &
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
 
   end function gmres_device_solve
 

src/krylov/bcknd/device/pipecg_device.F90

@@ -467,6 +467,7 @@ function pipecg_device_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(
       call this%monitor_stop()
       ksp_results%res_final = rnorm
       ksp_results%iter = iter
+      ksp_results%converged = this%is_converged(iter, rnorm)
 
     end associate
 

src/krylov/bcknd/sx/cg_sx.f90

@@ -34,7 +34,7 @@
 module cg_sx
   use num_types, only: rp
   use krylov, only : ksp_t, ksp_monitor_t, KSP_MAX_ITER
-  use precon,  only : pc_t
+  use precon, only : pc_t
   use ax_product, only : ax_t
   use field, only : field_t
   use coefs, only : coef_t
@@ -198,6 +198,7 @@ function sx_cg_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_resu
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
   end function sx_cg_solve
 
   !> Standard PCG coupled solve
@@ -220,9 +221,9 @@ function sx_cg_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
     type(ksp_monitor_t), dimension(3) :: ksp_results
     integer, optional, intent(in) :: niter
 
-    ksp_results(1) =  this%solve(Ax, x, fx, n, coef, blstx, gs_h, niter)
-    ksp_results(2) =  this%solve(Ax, y, fy, n, coef, blsty, gs_h, niter)
-    ksp_results(3) =  this%solve(Ax, z, fz, n, coef, blstz, gs_h, niter)
+    ksp_results(1) = this%solve(Ax, x, fx, n, coef, blstx, gs_h, niter)
+    ksp_results(2) = this%solve(Ax, y, fy, n, coef, blsty, gs_h, niter)
+    ksp_results(3) = this%solve(Ax, z, fz, n, coef, blstz, gs_h, niter)
 
   end function sx_cg_solve_coupled
 

src/krylov/bcknd/sx/gmres_sx.f90

@@ -33,7 +33,7 @@
 !> Defines various GMRES methods
 module gmres_sx
   use krylov, only : ksp_t, ksp_monitor_t
-  use precon,  only : pc_t
+  use precon, only : pc_t
   use ax_product, only : ax_t
   use num_types, only: rp
   use field, only : field_t
@@ -71,7 +71,7 @@ module gmres_sx
 
   !> Initialise a standard GMRES solver
   subroutine sx_gmres_init(this, n, max_iter, M, lgmres, &
-                           rel_tol, abs_tol, monitor)
+       rel_tol, abs_tol, monitor)
     class(sx_gmres_t), intent(inout) :: this
     integer, intent(in) :: n
     integer, intent(in) :: max_iter
@@ -200,7 +200,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
     integer :: i, j, k, ierr
     real(kind=rp), parameter :: one = 1.0
     real(kind=rp) :: rnorm
-    real(kind=rp) ::  alpha, temp, l
+    real(kind=rp) :: alpha, temp, l
     real(kind=rp) :: ratio, div0, norm_fac
     logical :: conv
     integer outer
@@ -232,7 +232,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
           call col3(this%r,this%ml,f,n)
        else
           !update residual
-          call copy  (this%r,f,n)
+          call copy (this%r,f,n)
           call Ax%compute(this%w, x%x, coef, x%msh, x%Xh)
           call gs_h%op(this%w, n, GS_OP_ADD)
           call bc_list_apply(blst, this%w, n)
@@ -284,7 +284,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
           !apply Givens rotations to new column
           do i=1,j-1
              temp = this%h(i,j)
-             this%h(i  ,j) =  this%c(i)*temp + this%s(i)*this%h(i+1,j)
+             this%h(i ,j) = this%c(i)*temp + this%s(i)*this%h(i+1,j)
              this%h(i+1,j) = -this%s(i)*temp + this%c(i)*this%h(i+1,j)
           end do
 
@@ -297,10 +297,10 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
           l = sqrt(this%h(j,j) * this%h(j,j) + alpha**2)
           temp = one / l
           this%c(j) = this%h(j,j) * temp
-          this%s(j) = alpha  * temp
+          this%s(j) = alpha * temp
           this%h(j,j) = l
           this%gam(j+1) = -this%s(j) * this%gam(j)
-          this%gam(j)   =  this%c(j) * this%gam(j)
+          this%gam(j) = this%c(j) * this%gam(j)
 
           rnorm = abs(this%gam(j+1)) * norm_fac
           call this%monitor_iter(iter, rnorm)
@@ -336,6 +336,7 @@ function sx_gmres_solve(this, Ax, x, f, n, coef, blst, gs_h, niter) result(ksp_r
     call this%monitor_stop()
     ksp_results%res_final = rnorm
     ksp_results%iter = iter
+    ksp_results%converged = this%is_converged(iter, rnorm)
   end function sx_gmres_solve
 
   !> Standard GMRES coupled solve
@@ -358,9 +359,9 @@ function sx_gmres_solve_coupled(this, Ax, x, y, z, fx, fy, fz, &
     type(ksp_monitor_t), dimension(3) :: ksp_results
     integer, optional, intent(in) :: niter
 
-    ksp_results(1) =  this%solve(Ax, x, fx, n, coef, blstx, gs_h, niter)
-    ksp_results(2) =  this%solve(Ax, y, fy, n, coef, blsty, gs_h, niter)
-    ksp_results(3) =  this%solve(Ax, z, fz, n, coef, blstz, gs_h, niter)
+    ksp_results(1) = this%solve(Ax, x, fx, n, coef, blstx, gs_h, niter)
+    ksp_results(2) = this%solve(Ax, y, fy, n, coef, blsty, gs_h, niter)
+    ksp_results(3) = this%solve(Ax, z, fz, n, coef, blstz, gs_h, niter)
 
   end function sx_gmres_solve_coupled
 

src/krylov/krylov.f90

@@ -60,6 +60,8 @@ module krylov
      real(kind=rp) :: res_start
      !> FInal residual
      real(kind=rp) :: res_final
+     !> Status
+     logical :: converged = .false.
   end type ksp_monitor_t
 
   !> Base abstract type for a canonical Krylov method, solving \f$ Ax = f \f$.
@@ -87,6 +89,8 @@ module krylov
      procedure, pass(this) :: monitor_stop => krylov_monitor_stop
      !> Monitor iteration
      procedure, pass(this) :: monitor_iter => krylov_monitor_iter
+     !> Check for convergence
+     procedure, pass(this) :: is_converged => krylov_is_converged
      !> Destructor.
      procedure(ksp_t_free), pass(this), deferred :: free
   end type ksp_t
@@ -328,4 +332,24 @@ subroutine krylov_monitor_iter(this, iter, rnorm)
 
   end subroutine krylov_monitor_iter
 
+  !> Check for convergence
+  !!
+  !! This function checks if the Krylov solver has converged.
+  !! The solver is considered converged if the residual is less than the
+  !! absolute tolerance.
+  !!
+  !! @param residual Residual
+  !! @param iter Iteration number
+  pure function krylov_is_converged(this, iter, residual) result(converged)
+    class(ksp_t), intent(in) :: this
+    integer, intent(in) :: iter
+    real(kind=rp), intent(in) :: residual
+    logical :: converged
+
+    converged = .true.
+    if (iter .ge. this%max_iter) converged = .false.
+    if (residual .gt. this%abs_tol) converged = .false.
+
+  end function krylov_is_converged
+
 end module krylov