Fixed petsc

Author: HarryMclean

src/heatflow/mod_evolve.f90

@@ -37,6 +37,9 @@ module evolution
 
   private
   public :: simulate
+  
+  ! Module-level variables for PETSc (persist across time steps)
+  integer, allocatable, save :: ia32(:), ja32(:)   ! 32-bit copies for PETSc
 
 contains
 
@@ -52,7 +55,6 @@ subroutine simulate(itime)
     integer:: ncg, itol, itmax !, iss
     integer :: iter
     real(real12) :: e, err, tol
-    integer, allocatable :: ia32(:), ja32(:)   ! 32-bit copies for PETSc
     integer :: NA32
 
     !----------------------
@@ -117,6 +119,16 @@ subroutine simulate(itime)
     !---------------------------------------------
     if ( iSteady .eq. 0 ) then
        S = - inverse_time * Temp_p * lin_rhoc - Qdens - B
+       if (IVERB .gt. 3) then
+          write(*,*) "S construction diagnostics:"
+          write(*,*) "  inverse_time =", inverse_time
+          write(*,*) "  Temp_p avg =", sum(Temp_p)/size(Temp_p)
+          write(*,*) "  lin_rhoc avg =", sum(lin_rhoc)/size(lin_rhoc)
+          write(*,*) "  Qdens avg =", sum(Qdens)/size(Qdens)
+          write(*,*) "  B avg =", sum(B)/size(B)
+          write(*,*) "  -inverse_time*Temp_p*lin_rhoc avg =", sum(-inverse_time*Temp_p*lin_rhoc)/size(Temp_p)
+          write(*,*) "  S before S_CAT avg =", sum(S)/size(S)
+       end if
        if ( iCAttaneo  .eq. 1) then
           S = S + S_CAT
        end if
@@ -159,6 +171,17 @@ subroutine simulate(itime)
    x = Temp_p + (Temp_p - Temp_pp)
    if (any(x - Temp_p .lt. TINY)) x = x + TINY ! avoid nan solver issue
 
+   ! Debug: Print initial guess statistics
+   if (IVERB .gt. 3) then
+      write(*,*) "========== PETSc Solver Diagnostics =========="
+      write(*,*) "Time step:", itime
+      write(*,*) "Initial guess x: min=", minval(x), " max=", maxval(x), " avg=", sum(x)/size(x)
+      write(*,*) "RHS S: min=", minval(S), " max=", maxval(S), " avg=", sum(S)/size(S)
+      write(*,*) "Temp_p: min=", minval(Temp_p), " max=", maxval(Temp_p), " avg=", sum(Temp_p)/size(Temp_p)
+      write(*,*) "Matrix acsr: min=", minval(acsr), " max=", maxval(acsr), " avg=", sum(acsr)/size(acsr)
+      write(*,*) "Matrix size: n=", NA32, " nnz=", size(acsr)
+   end if
+   
    ! Convert to 32-bit integers for PETSc (only on first call)
    if (.not. allocated(ia32)) then
       allocate(ia32(size(ia)), ja32(size(ja)))
@@ -169,6 +192,14 @@ subroutine simulate(itime)
 
    call solve_petsc_csr(NA32, ia32, ja32, acsr, S, x, tol, itmax)
 
+   ! Debug: Print solution statistics and verify solution
+   if (IVERB .gt. 3) then
+      write(*,*) "Solution x after PETSc: min=", minval(x), " max=", maxval(x), " avg=", sum(x)/size(x)
+      write(*,*) "Temperature change: avg(x-Temp_p)=", sum(x-Temp_p)/size(x)
+      write(*,*) "Max temperature change: ", maxval(abs(x-Temp_p))
+      write(*,*) "=============================================="
+   end if
+   
    ! Note: Don't deallocate ia32, ja32 - keep them for next time step
 
 

src/heatflow/mod_heating.f90

@@ -181,7 +181,7 @@ subroutine heater(itime, Q, Qdens)
 
 
     ! Normalize all heat sources by the heated volume
-    shared_power = .False.
+    shared_power = .TRUE.
     if (shared_power) then
       if (heated_volume .gt. 0.0) then
          Qdens(:) = Q(:) / heated_volume

src/heatflow/mod_petsc_solver.f90

@@ -42,19 +42,34 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
     ! Create matrix with an estimated 7 nonzeros/row (adjust if needed)
     call MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 7, PETSC_NULL_INTEGER, A, ierr)
 
+    ! Fill matrix from CSR format (ia, ja are 1-based Fortran indexing)
     do i = 1, n
        row_nz = ia(i+1) - ia(i)
        if (row_nz > 0) then
           start_k = ia(i)
           allocate(cols0(row_nz), vals(row_nz))
-          cols0 = ja(start_k:start_k+row_nz-1) - 1      ! zero-based
+          ! Convert column indices from 1-based to 0-based for PETSc
+          cols0 = ja(start_k:start_k+row_nz-1) - 1
           vals  = aval(start_k:start_k+row_nz-1)
+          
+          ! Debug: print first row details
+          if (i == 1) then
+             write(*,'(A,I0,A,I0)') 'Row 1: nnz=', row_nz, ', start_k=', start_k
+             write(*,'(A,10I6)') '  1-based cols:', ja(start_k:min(start_k+9,start_k+row_nz-1))
+             write(*,'(A,10I6)') '  0-based cols:', cols0(1:min(10,row_nz))
+             write(*,'(A,10ES12.4)') '  vals:', vals(1:min(10,row_nz))
+          end if
+          
+          ! Set row i-1 (0-based) with column indices cols0 (0-based)
           call MatSetValues(A, 1, (/i-1/), row_nz, cols0, vals, INSERT_VALUES, ierr)
           deallocate(cols0, vals)
        end if
     end do
     call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
     call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
+    
+    ! Optional: Verify matrix assembly (uncomment for debugging)
+    ! call MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr)
 
     ! Create vectors
     call VecCreateSeq(PETSC_COMM_SELF, n, bb, ierr)
@@ -93,6 +108,10 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
     call KSPGetIterationNumber(ksp, its, ierr)
     call KSPGetResidualNorm(ksp, rnorm, ierr)
 
+    ! Report convergence status (commented out by default for performance)
+    ! Uncomment the next line to see convergence info every solve:
+    ! write(*,'(A,I0,A,ES12.5,A,I0)') ' PETSc: iterations=', its, ', residual=', rnorm, ', reason=', reason
+    
     if (reason < 0) then
        write(0,*) "WARNING: PETSc solver diverged or failed!"
        write(0,*) "  Reason code:", reason
@@ -105,6 +124,37 @@ subroutine solve_petsc_csr(n, ia, ja, aval, b, x, rtol, maxit)
     call VecGetArrayF90(xx, xptr, ierr)
     x = xptr
     call VecRestoreArrayF90(xx, xptr, ierr)
+    
+    ! Manual residual check: compute ||A*x - b||
+    block
+      real(8), allocatable :: Ax(:), residual(:)
+      real(8) :: manual_rnorm, bnorm
+      integer(8) :: i, k, start_k, row_nz
+      
+      allocate(Ax(n), residual(n))
+      Ax = 0.0_8
+      
+      ! Compute A*x manually using CSR format
+      do i = 1, n
+         row_nz = ia(i+1) - ia(i)
+         if (row_nz > 0) then
+            start_k = ia(i)
+            do k = start_k, start_k + row_nz - 1
+               Ax(i) = Ax(i) + aval(k) * x(ja(k))
+            end do
+         end if
+      end do
+      
+      residual = Ax - b
+      manual_rnorm = sqrt(sum(residual**2))
+      bnorm = sqrt(sum(b**2))
+      
+      write(*,'(A,ES15.7)') ' Manual residual check: ||A*x-b|| = ', manual_rnorm
+      write(*,'(A,ES15.7)') '                        ||b||     = ', bnorm
+      write(*,'(A,ES15.7)') '                        ||A*x-b||/||b|| = ', manual_rnorm/bnorm
+      
+      deallocate(Ax, residual)
+    end block
 
     ! Cleanup
     deallocate(idx)

src/heatflow/mod_setup.f90

@@ -37,6 +37,7 @@ subroutine set_global_variables()
       integer(int12) :: ix,iy,iz,index
       real(real12) :: kappa,kappa3D,h_conv,heat_capacity,rho,sound_speed,tau, em
       real(real12), dimension(3) :: vel
+
       allocate(Temp_cur(nx, ny, nz))
       allocate(Temp_p(NA))
       allocate(Temp_pp(NA))
@@ -82,7 +83,6 @@ subroutine set_global_variables()
          ! Allocate the arrays to hold the H matrix in CSR format
          allocate(acsr(ra%len), ja(ra%len), ia(ra%n+1))
          CALL coo2csr(ra%n, ra%len, ra%val, ra%irow, ra%jcol, acsr, ja, ia)
-         ! print*, ra%val               
       end if
       !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -106,7 +106,7 @@ subroutine sparse_Hmatrix()
       if (Periodicz) len = len + 2*nx*ny
       ra%n = NA ! The number of rows in the H matrix
       ra%len = len ! The number of non-zero elements in the H matrix
-      ! Allocate the arrays to hold the H matrix in sparse row storage
+      ! Allocate the arrays to hold the H matrix in sparse storage
       allocate(ra%val(len), ra%irow(len), ra%jcol(len))
       ra%val(:)=0
       ra%irow(:)=-2
@@ -147,7 +147,6 @@ subroutine sparse_Hmatrix()
              ra%irow(count) = i ! The row of the H matrix
              ra%jcol(count) = j ! The column of the H matrix
              count = count + 1 ! The number of non-zero elements in the H matrix
-             H0=hmatrixfunc(j,i) ! The value of the H matrix
              ra%val(count) = H0 ! The value of the H matrix
              ra%irow(count) = j ! The row of the H matrix
              ra%jcol(count) = i ! The column of the H matrix

src/heatflow/mod_sparse_solver.f90

@@ -60,10 +60,14 @@ subroutine coo2csr( nrow, &
     integer(kind=8), dimension(nrow+1), intent(out) :: ia
 
     ! Local variables.
-    integer (kind=8) :: i, iad, j, k, k0
+    integer (kind=8) :: i, iad, j, k, k0, row_start
+    integer (kind=8) :: dup_count
     real(kind=8) :: x
+    logical :: found_dup
+    integer(kind=8), dimension(nrow+1) :: ia_save  ! Save original row starts
 
     ia(1:nrow+1) = 0
+    dup_count = 0
 
     ! determine the row lengths.
 
@@ -83,20 +87,45 @@ subroutine coo2csr( nrow, &
 
     end do
 
+    ! Save the original row pointers
+    ia_save = ia
+    
     ! go through the structure once more. fill in output matrix.
+    ! This version handles duplicate (i,j) entries by summing them.
 
     do k = 1, nnz
 
        i = ir(k)
        j = jc(k)
        x = a(k)
-       iad = ia(i)
-       acsr(iad) = x
-       ja(iad) = j
-       ia(i) = iad + 1
+       
+       ! Search for existing entry in this row with same column
+       found_dup = .false.
+       row_start = ia_save(i)
+       do iad = row_start, ia(i)-1
+          if (ja(iad) == j) then
+             ! Found duplicate - sum the values
+             acsr(iad) = acsr(iad) + x
+             found_dup = .true.
+             dup_count = dup_count + 1
+             exit
+          end if
+       end do
+       
+       if (.not. found_dup) then
+          ! New entry - insert it
+          iad = ia(i)
+          acsr(iad) = x
+          ja(iad) = j
+          ia(i) = iad + 1
+       end if
 
     end do
 
+    if (dup_count > 0) then
+       write(*,'(A,I0,A)') 'WARNING: COO->CSR found and summed ', dup_count, ' duplicate entries'
+    end if
+    
     ! shift back ia.
 
     do j = nrow, 1, -1