Working petsc, not fully reproducing old yet

Author: HarryMclean

src/.vscode/settings.json

@@ -0,0 +1,5 @@
+{
+    "python-envs.defaultEnvManager": "ms-python.python:conda",
+    "python-envs.defaultPackageManager": "ms-python.python:conda",
+    "python-envs.pythonProjects": []
+}

src/heatflow/mod_evolve.f90

@@ -153,22 +153,30 @@ subroutine simulate(itime)
     err=E
 
    !  call bicgstab(acsr, ia, ja, S, itmax, Temp_p, x, iter)
-  if (.not. allocated(ia32)) then
-     allocate(ia32(size(ia)), ja32(size(ja)))
-     ia32 = int(ia, kind=kind(ia32))
-     ja32 = int(ja, kind=kind(ja32))
-  end if
-   NA32 = int(NA, kind=kind(NA32))
+   
+   ! Allocate and initialize x with a good initial guess
    allocate(x(NA))
+   x = Temp_p + (Temp_p - Temp_pp)
+   if (any(x - Temp_p .lt. TINY)) x = x + TINY ! avoid nan solver issue
+   
+   ! Convert to 32-bit integers for PETSc (only on first call)
+   if (.not. allocated(ia32)) then
+      allocate(ia32(size(ia)), ja32(size(ja)))
+      ia32 = int(ia, kind=kind(ia32))
+      ja32 = int(ja, kind=kind(ja32))
+   end if
+   NA32 = int(NA, kind=kind(NA32))
+   
    call solve_petsc_csr(NA32, ia32, ja32, acsr, S, x, tol, itmax)
-  if (allocated(ia32)) deallocate(ia32, ja32)
+   
+   ! Note: Don't deallocate ia32, ja32 - keep them for next time step
 
 
    ! CALL solve_pardiso(acsr, S, ia, ja, x)
    !  CALL linbcg(S,x,itol=int(itol,I4B),tol=tol, itmax=int(itmax,I4B), iter=iter, &
          ! err=E)
 
-   ! 
+   !
     if (any(isnan(x(:)))) then
        write(0,*) "fatal error: NAN in x tempurature vector"
        write(0,*) 'time step ', itime, "      T   ", sum(Temp_p)/size(Temp_p), E ,iter

src/heatflow/mod_petsc_solver.f90

@@ -29,10 +29,12 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
     Vec :: bb, xx
     KSP :: ksp
     PC  :: pc
-    integer :: ierr, i, row_nz, start_k
+    integer :: ierr, i, row_nz, start_k, its
     integer, allocatable :: cols0(:), idx(:)
     real(8), allocatable :: vals(:)
     real(8), pointer :: xptr(:)
+    KSPConvergedReason :: reason
+    real(8) :: rnorm
 
     if (size(ia) /= n+1) stop 'solve_petsc_csr: ia size mismatch'
     if (size(b) /= n .or. size(x) /= n) stop 'solve_petsc_csr: vector size mismatch'
@@ -71,12 +73,33 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
     call KSPCreate(PETSC_COMM_SELF, ksp, ierr)
     call KSPSetOperators(ksp, A, A, ierr)     ! 3-arg form (reuse automatically)
     call KSPGetPC(ksp, pc, ierr)
-    call PCSetType(pc, PCILU, ierr)           ! Override at runtime with -pc_type
-    call KSPSetType(ksp, KSPCG, ierr)         ! Use -ksp_type bcgs if not SPD
+    call PCSetType(pc, PCJACOBI, ierr)        ! Use Jacobi (diagonal) preconditioner to match linbcg
+    call KSPSetType(ksp, KSPBCGS, ierr)       ! Use BCGS to match linbcg behavior
+    
+    ! Set convergence tolerances
+    ! rtol = relative tolerance, atol = absolute tolerance (use default), dtol = divergence tolerance, maxits = max iterations
     call KSPSetTolerances(ksp, rtol, PETSC_DEFAULT_REAL, PETSC_DEFAULT_REAL, maxit, ierr)
+    
+    ! Use unpreconditioned norm (matching linbcg with itol=1)
+    call KSPSetNormType(ksp, KSP_NORM_UNPRECONDITIONED, ierr)
+    
+    ! Allow command line override of solver options
     call KSPSetFromOptions(ksp, ierr)
 
     call KSPSolve(ksp, bb, xx, ierr)
+    
+    ! Check convergence
+    call KSPGetConvergedReason(ksp, reason, ierr)
+    call KSPGetIterationNumber(ksp, its, ierr)
+    call KSPGetResidualNorm(ksp, rnorm, ierr)
+    
+    if (reason < 0) then
+       write(0,*) "WARNING: PETSc solver diverged or failed!"
+       write(0,*) "  Reason code:", reason
+       write(0,*) "  Iterations:", its
+       write(0,*) "  Residual norm:", rnorm
+       ! Don't stop - let the main code detect NaNs if needed
+    end if
 
     ! Extract solution
     call VecGetArrayF90(xx, xptr, ierr)