Add ILU, and fix heating case 2

Author: HarryMclean

src/heatflow.f90

@@ -34,6 +34,12 @@ program HEATFLOW_V0_3
    real(real12) :: cpustart, cpuend, cpustart2, progress
    integer(int12) :: itime
 
+   !-------------------------------------------------------------!
+   ! Initialize PETSc FIRST (before any other operations)       !
+   !-------------------------------------------------------------!
+   CALL petsc_init()
+   !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
+
    !-------------------------------------------------------------!
    ! calculate the time to run full simulation                   !
    !-------------------------------------------------------------!
@@ -65,7 +71,6 @@ program HEATFLOW_V0_3
    !-------------------------------------------------------------!
    ! run simulation for 'ntime' time steps                       !
    !-------------------------------------------------------------!
-   CALL petsc_init()                                 
 
    do itime=1,ntime 
 
@@ -99,7 +104,7 @@ program HEATFLOW_V0_3
    !-------------------------------------------------------------!
    ! calculate end time and print to user                        !
    !-------------------------------------------------------------!
-   CALL cpu_time(cpuend)                                          
+   CALL cpu_time(cpuend)
    write(*,'(A,F12.6)') ' time=', cpuend-cpustart                 
    !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
 

src/heatflow/mkl_pardiso.f90

@@ -159,7 +159,7 @@ subroutine simulate(itime)
    !  if (any(x-Temp_p .lt. TINY)) x=x+TINY !avoid nan solver issue
     itol=1
     tol=1.e-32_real12
-    itmax=50000
+    itmax=500000
     ncg = 0
     iter= 0
     err=E

src/heatflow/mod_evolve.f90

@@ -78,18 +78,22 @@ subroutine heater(itime, Q, Qdens)
                 !------------------------------
                 if ( time .le. time_pulse ) then
                    Q(IA) = POWER
-                if (icattaneo .eq. 1) then
-                if (itime .eq. 1) then
-                  Q(IA) = Q(IA) + (tau*(POWER))
-                  end if
-                if (itime .eq. heated_steps+1) then
-                  Q(IA) = Q(IA) - (tau*(POWER))
-                end if
-                end if 
-                
+                  !  print *, "Heating on"
                 else
+                  !  print *, "Heating off"
                    Q(IA) = 0.0_real12
                 end if
+
+               if (icattaneo .eq. 1) then
+                  if (itime .eq. 1) then
+                     Q(IA) = Q(IA) + tau*POWER
+                     ! print *, "Turning on heater at time ", time
+                  end if
+                  if (itime .eq. (heated_steps + 1)) then
+                     Q(IA) = Q(IA) - tau*POWER
+                     ! print *, "Turning off heater at time ", time
+                  end if
+               end if
 
 
                 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -126,6 +130,34 @@ subroutine heater(itime, Q, Qdens)
                 end if
                 !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
+            case(7)
+               !---------------------------------------------------------
+               ! Heater on for a time period, off for a time period,
+               ! then on again (square-wave heating)
+               !---------------------------------------------------------
+
+               if (mod(time, 2.0_real12 * time_pulse) .le. time_pulse) then
+                  Q(IA) = POWER
+
+                  if (icattaneo .eq. 1) then
+
+                     ! Forward-time pulse check
+                     if (mod(time + time_step, 2.0_real12 * time_pulse) .le. time_pulse) then
+                        Q(IA) = Q(IA) + tau * POWER
+                     end if
+
+                     ! Backward-time pulse check
+                     if (mod(time - time_step, 2.0_real12 * time_pulse) .le. time_pulse) then
+                        Q(IA) = Q(IA) - tau * POWER
+                     end if
+
+                  end if
+
+               else
+                  Q(IA) = 0.0_real12
+               end if
+               
+
                case(10)
                   Q(IA) = POWER - (tau*(POWER))