diff --git a/advanced_lab/c/shwater2d.c b/advanced_lab/c/shwater2d.c index 0d87f665fe5ea89838705d6689fab77f3c44aad3..90b905c1312b3d197487d08a0871be36d6e0d981 100644 --- a/advanced_lab/c/shwater2d.c +++ b/advanced_lab/c/shwater2d.c @@ -1,5 +1,5 @@ /* - * 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 */ diff --git a/advanced_lab/f90/shwater2d.f90 b/advanced_lab/f90/shwater2d.f90 index 5e651fcf1b3a58129e9a13ec86e0c013d4f725cc..9ec9d1fd63f0794d4957d7fb1bd327e5c1d33ea8 100644 --- a/advanced_lab/f90/shwater2d.f90 +++ b/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,8 +106,15 @@ 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 + +!-----------------------------------------------------------------------