Skip to content
Snippets Groups Projects
Commit 24a24061 authored by Kjartan Thor Wikfeldt's avatar Kjartan Thor Wikfeldt
Browse files

Merge branch 'master' of github.com:PDC-support/openmp-lab-exercises

parents b95d9fe5 774836da
No related branches found
No related tags found
No related merge requests found
/*
* 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 */
......
......@@ -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
!-----------------------------------------------------------------------
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment