Merge pull request #172 from BerkeleyLab/fix-examples

fix(examples): time step stability calculation

Author: rouson

example/heat-equation.f90

@@ -19,8 +19,7 @@ module subdomain_2D_m
     real, allocatable :: s_(:,:)
   contains
     procedure define
-    procedure dx
-    procedure dy
+    procedure dt_stable
     procedure values
     procedure exchange_halo
     generic :: operator(.laplacian.) => laplacian
@@ -44,16 +43,11 @@ pure module function values(self) result(my_values)
       real, allocatable :: my_values(:,:)
     end function
 
-    pure module function dx(self) result(my_dx)
+    pure module function dt_stable(self, alpha) result(my_dt)
       implicit none
       class(subdomain_2D_t), intent(in) :: self
-      real my_dx
-    end function
-
-    pure module function dy(self) result(my_dy)
-      implicit none
-      class(subdomain_2D_t), intent(in) :: self
-      real my_dy
+      real, intent(in) :: alpha
+      real my_dt
     end function
 
     pure module function laplacian(rhs) result(laplacian_rhs)
@@ -118,12 +112,11 @@ module subroutine exchange_halo(self)
     call self%exchange_halo
   end procedure
 
-  module procedure dx
-    my_dx = dx_
-  end procedure
-
-  module procedure dy
-    my_dy = dy_
+  module procedure dt_stable
+    !! Set the time step at 90% of the stability limit obtained provided by
+    !! Kassinos, S., & Alexiadis, A. (2024). Beyond Language: Applying MLX Transformers to Engineering Physics
+    !! arXiv preprint arXiv:2410.04167.
+    my_dt = 0.9 * (1./(1./dx_**2 + 1./dy_**2)) *  (1./(2.*alpha))
   end procedure
 
   module procedure laplacian
@@ -192,7 +185,7 @@ program main
   call T%define(side=1., boundary_val=T_boundary, internal_val=T_initial, n=nx) ! 2D step function
   sync all
 
-  associate(dt => T%dx()*T%dy()/(4*alpha))
+  associate(dt => T%dt_stable(alpha))
     call cpu_time(t_start)
     do step = 1, steps
       T =  T + dt * alpha * .laplacian. T

example/time-paradigm.f90

@@ -58,7 +58,7 @@ function functional_programming_time() result(system_time)
 
     call system_clock(t_start_functional)
 
-    associate(dt => T%dx()*T%dy()/(4*alpha))
+    associate(dt => T%dt_stable(alpha))
       functional_programming: &
       do step = 1, steps
         T =  T + dt * alpha * .laplacian. T
@@ -80,7 +80,7 @@ function procedural_programming_time() result(system_time)
     real system_time
     type(subdomain_t) T
 
-    associate(dt => T%dx()*T%dy()/(4*alpha))
+    associate(dt => T%dt_stable(alpha))
       call T%define(side=1., boundary_val=0., internal_val=1., n=resolution)
       call system_clock(t_start_procedural)
       procedural_programming: &

test/subdomain_test_m.f90

@@ -160,7 +160,7 @@ function correct_steady_state() result(test_passes)
 
     call T%define(side=1., boundary_val=T_boundary, internal_val=T_initial, n=11) ! const. internally with a step down at boundaries
 
-    associate(dt => t%dt_stable(alpha)) 
+    associate(dt => T%dt_stable(alpha))
       do step = 1, steps
         T =  T + dt * alpha * .laplacian. T
       end do
@@ -193,7 +193,7 @@ function T_functional()
 
       call T%define(side, boundary_val, internal_val, n)
 
-      associate(dt => t%dt_stable(alpha)) 
+      associate(dt => T%dt_stable(alpha))
         do step = 1, steps
           T =  T + dt * alpha * .laplacian. T
         end do
@@ -209,7 +209,7 @@ function T_procedural()
 
       call T%define(side, boundary_val, internal_val, n)
 
-      associate(dt => t%dt_stable(alpha)) 
+      associate(dt => T%dt_stable(alpha))
         do step = 1, steps
           call T%step(alpha*dt)
         end do