Add rudimentary error checks

advanced_lab/c/shwater2d.c

 /*
- *  shwater2d.f90 solves the two dimensional shallow water equations 
+ *  shwater2d.c solves the two dimensional shallow water equations 
  *  using the Lax-Friedrich's scheme
  */
 
@@ -24,6 +24,18 @@ double gettime(void) {
   return tv.tv_sec + 1e-6*tv.tv_usec;
 }
 
+/* Check that the solution is finite */
+void validate(double *Q, int m, int n) {
+  int i, j, k;
+  for (i = 0; i < n; i++) 
+    for (j = 0; j < m; j++) 
+      for (k = 0; k < cell_size;  k++)
+	if (!isfinite(Q(k, j, i))) {
+	  fprintf(stderr, "Invalid solution\n");
+	  exit(-1);
+	}
+}
+
 /* Flux function in the x-direction */
 void fx(double *Q, double **fq, int m, int n, int j) {  
   int i;
@@ -136,7 +148,7 @@ int main(int argc, char **argv) {
   double *Q;
   double *x, *y;
   double **ffx, **nFx, **ffy, **nFy;
-  double dx, dt, epsi, delta, dy, tend, tmp, stime;
+  double dx, dt, epsi, delta, dy, tend, tmp, stime, etime;
   
 
   /* Use m volumes in the x-direction and n volumes in the y-direction */    
@@ -210,7 +222,11 @@ int main(int argc, char **argv) {
 
   stime = gettime();
   solver(Q, ffx, ffy, nFx, nFy, m, n, tend, dx, dy, dt);
-  printf("Solver took %g seconds\n", gettime() - stime);
+  etime = gettime();
+
+  validate(Q, m,  n);
+
+  printf("Solver took %g seconds\n", etime - stime);
 
 
   /* Uncomment this line if you want visualize the result in ParaView */

advanced_lab/f90/shwater2d.f90

@@ -35,7 +35,7 @@ program shwater2d
   real(kind=dp), dimension(:), allocatable :: x, y
   integer i, j, ifail, m, n
   real(kind=dp) dx, dt, epsi, delta, dy, tend, tmp
-  real stime
+  real stime, etime
   real, external :: rtc
   external fx,fy
 
@@ -73,8 +73,7 @@ program shwater2d
 
   if(ifail .ne. 0) then
      deallocate(Q, x, y)
-     write(*,*) 'Memory exhausted'
-     stop
+     stop 'Memory exhausted'
   else
      tmp = -dx/2 + xstart
      do i=1,m
@@ -107,7 +106,14 @@ program shwater2d
   !
   stime = rtc()
   call solver(Q, cell_size, m, n, tend, dx, dy, dt, fx, fy)
-  write(*,*) 'Solver took',rtc()-stime,'sec'
+  etime = rtc()
+
+  !
+  ! Check if solution is finite
+  ! 
+  call validate(Q, cell_size, m, n)
+
+  write(*,*) 'Solver took',etime-stime,'sec'
   
   !
   ! Uncomment this line if you want visualize the result in ParaView
@@ -229,6 +235,7 @@ end subroutine laxf_scheme_2d
 ! fx              flux function in the x-direction
 ! fy              flux function in the y-direction
 ! rtc             timing function
+! validate        check that the solution is finite
 !
 !-----------------------------------------------------------------------
 
@@ -284,3 +291,25 @@ real function rtc()
 end function rtc
 
 !-----------------------------------------------------------------------
+
+subroutine validate(q, l, m, n)
+  use types
+  use, intrinsic :: ieee_arithmetic, only: ieee_is_finite
+  implicit none
+  real(kind=dp), intent(in),dimension(l,m,n) :: q
+  integer, intent(in) :: l, m, n
+  integer i, j, k 
+
+  do j=1,n
+     do i=1,m
+        do k=1,l
+           if (.not. ieee_is_finite(q(k,i,j))) then  
+              stop 'Invalid solution'
+           end if
+        end do
+     end do
+  end do
+
+end subroutine validate
+
+!-----------------------------------------------------------------------