add advanced lab from unknown year

advanced_lab/c/Makefile

+CC = cc
+CFLAGS = -O2 -mp
+SRC = vtk_export.c shwater2d.c
+OBJS = ${SRC:.c=.o}
+DEST = shwater2d
+
+all: $(DEST)
+
+$(DEST): $(OBJS)
+	$(CC) $(CFLAGS) $(OBJS) -o $@ -lm	
+
+clean:
+	rm -f $(DEST) *.mod *.MOD *.o
+
+%.o: %.c
+	$(CC) $(CFLAGS) -c $<
+

advanced_lab/c/shwater2d.c

+/*
+ *  shwater2d.f90 solves the two dimensional shallow water equations 
+ *  using the Lax-Friedrich's scheme
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+//#include <omp.h>
+#include <sys/time.h>
+
+#define cell_size 3
+#define xstart 0.0
+#define ystart 0.0
+#define xend 4.0
+#define yend 4.0
+
+#define Q(i, j, k) Q[((k) + n * ((j) + m * (i)))]
+
+/* Timing function */
+double gettime(void) {
+  struct timeval tv;
+  gettimeofday(&tv,NULL);
+  return tv.tv_sec + 1e-6*tv.tv_usec;
+}
+
+/* Flux function in the x-direction */
+void fx(double *Q, double **fq, int m, int n, int j) {  
+  int i;
+  const double g = 9.81;
+
+  for (i = 0; i < m; i++) {
+    fq[0][i] = Q(1, i, j);
+    fq[1][i] = (pow(Q(1, i, j), 2) / Q(0, i, j))  + 
+      (g * pow(Q(0, i, j), 2)) / 2.0;
+    fq[2][i] = (Q(1, i, j) * Q(2, i, j)) / Q(0, i, j);
+  }
+  
+}
+
+/* Flux function in the y-direction */
+void fy(double *Q, double **fq, int m, int n, int i) {
+  int j;
+  const double g= 9.81;
+
+  for (j = 0; j < n; j++) {
+    fq[0][j] = Q(2, i, j);
+    fq[1][j] = (Q(1, i, j) * Q(2, i, j)) / Q(0, i, j);    
+    fq[2][j] = (pow(Q(2, i, j), 2) / Q(0, i, j))  + 
+      (g * pow(Q(0, i, j), 2)) / 2.0;
+  }
+}
+
+
+/*
+  This is the Lax-Friedrich's scheme for updating volumes
+  Try to parallelize it in an efficient way!
+*/
+void laxf_scheme_2d(double *Q, double **ffx, double **ffy, double **nFx, double **nFy,
+		    int m, int n, double dx, double dy, double dt) {
+  int i, j, k;
+    
+  /* Calculate and update fluxes in the x-direction */
+  for (i = 1; i < n; i++) {
+    fx(Q, ffx, m, n, i);
+    for (j = 1; j < m; j++) 
+      for (k = 0; k < cell_size;  k++) 
+	nFx[k][j] = 0.5 * ((ffx[k][j-1] + ffx[k][j]) -
+			   dx/dt * (Q(k, j, i) - Q(k, j-1, i)));
+    for (j = 1; j < m-1; j++)
+      for (k = 0; k < cell_size;  k++) 
+	Q(k, j, i) = Q(k, j, i)  - dt/dx * ((nFx[k][j+1] - nFx[k][j]));
+      
+  }
+
+  /* Calculate and update fluxes in the y-direction */
+  for (i = 1; i < m; i++) {
+    fy(Q, ffy, m, n, i);
+    for (j = 1; j < n; j++)
+      for (k = 0; k < cell_size; k++)
+	nFy[k][j] = 0.5 * ((ffy[k][j-1] + ffy[k][j]) - 
+			   dy/dt * (Q(k, i, j) - Q(k, i, j -1)));
+    for (j = 1; j <  n-1; j++) 
+      for (k = 0; k < cell_size; k++)
+	Q(k,i,j) = Q(k,i,j) -  dt/dy * ((nFy[k][j+1]  -  nFy[k][j]));
+  }
+
+}
+
+/*
+  This is the main solver routine, parallelize this. 
+  But don't forget the subroutine laxf_scheme_2d
+*/
+void solver(double *Q, double **ffx, double **ffy, double **nFx, double **nFy,
+	    int m, int n, double tend, double dx, double dy, double dt) {
+  double bc_mask[3] = {1.0, -1.0, -1.0};
+  double time;
+  int i, j, k, steps;
+  
+  steps = ceil(tend / dt);  
+  for (i = 0, time = 0.0; i < steps; i++, time += dt) { 
+
+    /* Apply boundary condition */
+    for (j = 1; j < n - 1 ; j++) {
+      for (k = 0; k < cell_size; k++) {
+	Q(k, 0, j) = bc_mask[k] *  Q(k, 1, j);
+	Q(k, m-1, j) = bc_mask[k] *  Q(k, m-2, j);
+      }
+    }
+
+    for (j = 0; j < m; j++)  {
+      for (k = 0; k < cell_size; k++) {
+	Q(k, j, 0) = bc_mask[k] * Q(k, j, 1);
+	Q(k, j, n-1) = bc_mask[k] * Q(k, j, n-2);
+      }
+    }
+     
+    /* Update all volumes with the Lax-Friedrich's scheme */     
+    laxf_scheme_2d(Q, ffx, ffy, nFx, nFy, m, n, dx, dy, dt);  
+  
+  }
+}
+  
+/*
+  This is the main routine of the program, which allocates memory 
+  and setup all parameters for the problem.
+  
+  You don't need to parallelize anything here!
+  
+  However, it might be useful to change the m and n parameters 
+  during debugging
+*/
+int main(int argc, char **argv) {
+  
+  int i, j, m, n;
+  double *Q;
+  double *x, *y;
+  double **ffx, **nFx, **ffy, **nFy;
+  double dx, dt, epsi, delta, dy, tend, tmp, stime;
+  
+
+  /* Use m volumes in the x-direction and n volumes in the y-direction */    
+  m = 1000;
+  n = 1000;
+  
+  /*
+    epsi      Parameter used for initial condition
+    delta     Parameter used for initial condition
+    dx        Distance between two volumes (x-direction)
+    dy        Distance between two volumes (y-direction)
+    dt        Time step
+    tend      End time
+  */
+  epsi = 2.0;
+  delta = 0.5;
+  dx = (xend - xstart) / (double) m;
+  dy = (yend - ystart) / (double) n;
+  dt = dx / sqrt( 9.81 * 5.0);
+  tend = 0.1;
+  
+  /* Add two ghost volumes at each side of the domain */
+  m = m + 2;
+  n = n + 2;
+
+  /* Allocate memory for the domain */
+  Q = (double *) malloc(m * n * cell_size *  sizeof(double));
+
+  x = (double *) malloc(m * sizeof(double));
+  y = (double *) malloc(n * sizeof(double));	
+
+  /* Allocate memory for fluxes */
+  ffx = (double **) malloc(cell_size * sizeof(double *));
+  ffy = (double **) malloc(cell_size * sizeof(double *));
+  nFx = (double **) malloc(cell_size * sizeof(double *));
+  nFy = (double **) malloc(cell_size * sizeof(double *));
+
+  ffx[0] = (double *) malloc(cell_size * m * sizeof(double));
+  ffy[0] = (double *) malloc(cell_size * m * sizeof(double));
+  nFx[0] = (double *) malloc(cell_size * n * sizeof(double));
+  nFy[0] = (double *) malloc(cell_size * n * sizeof(double));
+
+  for (i = 0; i < cell_size; i++) {
+    ffx[i] =  ffx[0] + i * m;
+    nFx[i] =  nFx[0] + i * m;
+    ffy[i] =  ffy[0] + i * n;
+    nFy[i] =  nFy[0] + i * n;
+  }
+    
+  for (i = 0,tmp= -dx/2 + xstart; i < m; i++, tmp += dx)
+    x[i] = tmp;
+
+  for (i = 0,tmp= -dy/2 + ystart; i < n; i++, tmp += dy)
+    y[i] = tmp;
+
+  /* Set initial Gauss hump */
+  for (i = 0; i < m; i++) {
+    for (j = 0; j < n; j++) {
+      Q(0, i, j) = 4.0;
+      Q(1, i, j) = 0.0;
+      Q(2, i, j) = 0.0;
+    }
+  }
+
+  for (i = 1; i < m-1; i++) {
+    for (j = 1; j < n-1; j++) {
+      Q(0, i, j) = 4.0 + epsi * exp(-(pow(x[i] - xend / 4.0, 2) + pow(y[j] - yend / 4.0, 2)) /
+					  (pow(delta, 2)));
+    }
+  }
+
+  stime = gettime();
+  solver(Q, ffx, ffy, nFx, nFy, m, n, tend, dx, dy, dt);
+  printf("Solver took %g seconds\n", gettime() - stime);
+
+
+  /* Uncomment this line if you want visualize the result in ParaView */
+  /* save_vtk(Q, x, y, m, n); */
+
+
+  free(Q);
+  free(x);
+  free(y);
+
+  free(ffx[0]);
+  free(ffy[0]);
+  free(nFx[0]);
+  free(nFy[0]);
+
+  free(ffx);
+  free(ffy);
+  free(nFx);
+  free(nFy);
+
+  return 0;    
+}
+
+
+

advanced_lab/c/vtk_export.c

+
+#include <stdlib.h>
+#include <stdio.h>
+
+#define Q(i, j, k) Q[((k) + (n) * ((j) + (m) * (i)))]
+
+
+
+void save_vtk(double *Q, double *x, double *y, int m, int n) {
+  
+  int i, j;
+  FILE *fp = fopen("result.vtk", "w");
+  
+  /* Write vtk Datafile header */
+  fprintf(fp, "# vtk DataFile Version 2.0\n");
+  fprintf(fp, "VTK\nASCII\nDATASET POLYDATA\n");
+
+  /* Store water height as polydata */
+  fprintf(fp, "\nPOINTS %d double\n", m*n);
+  
+  for (j = 0; j < n; j++) 
+    for (i = 0; i < m; i++)
+      fprintf(fp, "%e %e %e\n", x[i], y[j], Q(0, i, j));
+
+  fprintf(fp,"\nVERTICES %d %d\n", n, n *(m+1));
+  for (j = 0; j < n; j++)  {
+    fprintf(fp, "%d ", m);
+    for (i = 0; i < m; i++) 
+      fprintf(fp, "%d ", i+j*m);
+    fprintf(fp,"\n");
+  }
+
+  /* Store lookup table */
+  fprintf(fp,
+	  "POINT_DATA %d\nSCALARS height double 1\nLOOKUP_TABLE default\n",m*n);
+  for (j = 0; j < n; j++)
+    for (i = 0; i < m; i++)
+      fprintf(fp, "%e\n", Q(0, i, j));
+  fclose(fp);     
+}

advanced_lab/f90/Makefile

+F90 = ftn
+FFLAGS = -O2 -mp
+SRC = vtk_export.f90 shwater2d.f90
+OBJS = ${SRC:.f90=.o}
+DEST = shwater2d
+
+all: $(DEST)
+
+$(DEST): $(OBJS)
+	$(F90) $(FFLAGS) $(OBJS) -o $@	
+
+clean:
+	rm -f $(DEST) *.mod *.MOD *.o
+
+%.o: %.f90
+	$(F90) $(FFLAGS) -c $<
+

advanced_lab/f90/shwater2d.f90

+!-----------------------------------------------------------------------
+!
+! shwater2d.f90 solves the two dimensional shallow water equations 
+! using the Lax-Friedrich's scheme
+!
+!-----------------------------------------------------------------------
+
+module types
+  integer, parameter :: dp = kind(0.0d0)
+  integer, parameter :: sp = kind(0.0)
+end module types
+
+!-----------------------------------------------------------------------
+
+program shwater2d
+  use types
+  use omp_lib
+  use vtk_export
+  implicit none
+  
+  ! This is the main routine of the program, which allocates memory 
+  ! and setup all parameters for the problem.
+  !
+  ! You don't need to parallelize anything here!
+  !
+  ! However, it might be useful to change the m and n parameters 
+  ! during debugging
+  !
+   
+  integer, parameter :: cell_size = 3  
+  real(kind=dp) xstart, xend, ystart, yend
+  parameter (xstart = 0.0d0, ystart = 0.0d0, xend = 4d0, yend = 4d0)
+
+  real(kind=dp), dimension(:,:,:), allocatable :: Q
+  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, external :: rtc
+  external fx,fy
+
+
+  !
+  ! Use m volumes in the x-direction and n volumes in the y-direction
+  !
+  m = 1000
+  n = 1000
+  
+
+  ! epsi      Parameter used for initial condition
+  ! delta     Parameter used for initial condition
+  ! dx        Distance between two volumes (x-direction)
+  ! dy        Distance between two volumes (y-direction)
+  ! dt        Time step
+  ! tend      End time
+  epsi = 2d0
+  delta = 0.5
+  dx = (xend - xstart) / m
+  dy = (yend - ystart) / n
+  dt = dx / sqrt( 9.81d0 * 5d0) 
+  tend = 0.1
+
+  !
+  ! Add two ghost volumes at each side of the domain
+  !
+  m = m + 2
+  n = n + 2
+
+  !
+  ! Allocate memory for the domain
+  !
+  allocate(Q(cell_size, m, n), x(m), y(n), stat = ifail)
+
+  if(ifail .ne. 0) then
+     deallocate(Q, x, y)
+     write(*,*) 'Memory exhausted'
+     stop
+  else
+     tmp = -dx/2 + xstart
+     do i=1,m
+        x(i) = tmp
+        tmp = tmp + dx
+     end do
+
+     tmp = -dy/2 + ystart
+     do i=1,n
+        y(i) = tmp
+        tmp = tmp + dy
+     end do     
+  end if
+
+  !
+  ! Set initial Gauss hump
+  !
+  Q(2,:,:) = 0
+  Q(3,:,:) = 0
+  Q(1,:,:) = 4
+  do j=2,n-1
+     do i=2,m-1
+        Q(1,i,j) =  4 + epsi * exp(-((x(i) - xend / 4d0)**2 + &
+             (y(j) - yend / 4d0)**2)/delta**2)
+     end do
+  end do
+
+  !
+  ! Start solver
+  !
+  stime = rtc()
+  call solver(Q, cell_size, m, n, tend, dx, dy, dt, fx, fy)
+  write(*,*) 'Solver took',rtc()-stime,'sec'
+
+  !
+  ! Uncomment this line if you want visualize the result in ParaView
+  !
+  ! call save_vtk(Q, x, y, cell_size, m, n, 'result.vtk')
+
+  deallocate(Q, x, y)
+
+end program shwater2d
+
+!-----------------------------------------------------------------------
+
+subroutine solver(Q, l, m, n, tend, dx, dy, dt, fx, fy)
+  use types
+  implicit none
+
+  !
+  ! This is the main solver routine, parallelize this. 
+  ! But don't forget the subroutine laxf_scheme_2d
+  !
+
+  integer, intent(in) :: l, m, n
+  real(kind=dp), intent(inout), dimension(l, m, n) :: Q
+  real(kind=dp), intent(in) :: dx, dy, dt
+  real(kind=dp), dimension(l) :: bc_mask
+  real(kind=dp), intent(in) :: tend
+  real(kind=dp)  :: time
+  external fx, fy
+  integer i, j, steps
+
+  bc_mask(1) = 1d0
+  bc_mask(2:l) = -1d0
+
+  steps = tend / dt
+  time = 0d0
+
+  !
+  ! This is the main time stepping loop
+  !
+  do i=1,steps
+     
+     !
+     ! Apply boundary condition
+     !
+     do j  = 2, n - 1
+        Q(:, 1, j) = bc_mask * Q(:, 2, j)
+        Q(:, m, j) = bc_mask * Q(:, m-1, j)
+     end do
+
+     do j  = 1, m
+        Q(:, j, 1) = bc_mask *  Q(:, j , 2)
+        Q(:, j, n) = bc_mask *  Q(:, j, n-1)
+     end do
+     
+     !
+     ! Update all volumes with the Lax-Friedrich's scheme
+     !
+     call laxf_scheme_2d(Q, l, m, n, dx, dy, dt, fx, fy)
+     time = time + dt
+
+  end do
+
+end subroutine solver
+
+!-----------------------------------------------------------------------
+
+subroutine laxf_scheme_2d(Q, l, m, n, dx, dy, dt, fx, fy)
+  use types
+  implicit none 
+
+  !
+  ! This is the Lax-Friedrich's scheme for updating volumes
+  ! Try to parallelize it in an efficient way!
+  ! 
+
+  integer, intent(in) :: l, m, n
+  real(kind=dp), intent(inout), dimension(l, m, n) :: Q
+  real(kind=dp), dimension(l, m) :: ffx, nFx
+  real(kind=dp), dimension(l, n) :: ffy, nFy
+  real(kind=dp), intent(in) :: dx, dy, dt
+  external fx, fy
+  integer i,j
+
+  !
+  ! Calculate and update fluxes in the x-direction
+  !
+  do i=2,n
+     call fx(Q(:,:,i), ffx, l, m)
+     do j=2,m
+        nFx(:,j) = 0.5  * ((ffx(:,j-1) + ffx(:,j)) - &
+             dx/dt * (Q(:,j,i) - Q(:,j-1,i)))
+     end do     
+     do j=2,m-1
+        Q(:,j,i) = Q(:,j,i) -  dt/dx * ((nFx(:,j+1)  -  nFx(:,j)))
+     end do
+  end do
+
+  !
+  ! Calculate and update fluxes in the y-direction
+  !
+  do i=2,m
+     call fy(Q(:,i,:), ffy, l, n)
+     do j=2,n
+        nFy(:,j) = 0.5  * ((ffy(:,j-1) + ffy(:,j)) - &
+             dy/dt * (Q(:,i,j) - Q(:,i,j-1)))
+     end do
+     do j=2,n-1
+        Q(:,i,j) = Q(:,i,j) -  dt/dy * ((nFy(:,j+1)  -  nFy(:,j)))
+     end do     
+  end do
+
+end subroutine laxf_scheme_2d
+
+!-----------------------------------------------------------------------
+!
+! The rest of the file contains auxiliary functions, which you don't
+! need to parallelize.
+!
+! fx              flux function in the x-direction
+! fy              flux function in the y-direction
+! rtc             timing function
+!
+!-----------------------------------------------------------------------
+
+subroutine fx(q, fq, m, n)
+  use types 
+  implicit none 
+
+  real(kind=dp) ,intent(in),dimension(m,n) :: q
+  real(kind=dp) ,intent(out),dimension(m,n) :: fq
+  integer ,intent(in) :: m,n
+  real(kind=dp) ,parameter :: g = 9.81
+
+  fq(1,:) = q(2,:)
+  fq(2,:) = (q(2,:)**2 / q(1,:)) + (g * q(1,:)**2) / 2
+  fq(3,:) = (q(2,:) * q(3,:)) / q(1,:)
+  
+end subroutine fx
+
+!-----------------------------------------------------------------------
+
+subroutine fy(q, fq, m, n)
+  use types
+  implicit none 
+
+  real(kind=dp) ,intent(in),dimension(m,n) :: q
+  real(kind=dp) ,intent(out),dimension(m,n) :: fq
+  integer ,intent(in) :: m,n
+  real(kind=dp) ,parameter :: g = 9.81
+
+  fq(1,:) = q(3,:)
+  fq(2,:) = (q(2,:) * q(3,:)) / q(1,:) 
+  fq(3,:) = (q(3,:)**2 / q(1,:) ) + (g * q(1,:)**2) / 2
+  
+end subroutine fy
+
+!-----------------------------------------------------------------------
+
+real function rtc()
+  implicit none
+  integer:: icnt,irate
+  real, save:: scaling
+  logical, save:: scale = .true.
+
+  call system_clock(icnt,irate)
+
+  if(scale)then
+     scaling=1.0/real(irate)
+     scale=.false.
+  end if
+
+  rtc = icnt * scaling
+
+end function rtc
+
+!-----------------------------------------------------------------------

advanced_lab/f90/vtk_export.f90

+module vtk_export
+  implicit none
+
+  ! VTK export
+  !
+  ! This module stores the volumes Q first variable 
+  ! as VTK polydata
+  !
+  
+contains
+  subroutine save_vtk(Q, x, y, l, m, n, fname, title)
+    implicit none
+    integer, intent(in) :: l, m, n
+    double precision, intent(in), dimension(l, m, n) :: Q
+    double precision, intent(in), dimension(n) :: x
+    double precision, intent(in), dimension(m) :: y
+    character, intent(in) :: fname*(*)
+    character, optional :: title*(*)
+    integer i,j
+
+    open(1, file=fname)
+
+    !
+    ! Write vtk Datafile header
+    !
+    write(1,fmt='(A)') '# vtk DataFile Version 2.0'
+     if(present(title)) then
+       write(1,fmt='(A)') title
+    else
+       write(1,fmt='(A)') 'VTK'
+    end if  
+    write(1,fmt='(A)') 'ASCII'
+    write(1,fmt='(A)') 'DATASET POLYDATA'
+    write(1,*)
+    
+    !
+    ! Store water height as polydata
+    !
+     write(1,fmt='(A,I8,A)') 'POINTS', m*n,' double'
+     do j=1,n
+        do i=1,m
+           write(1,fmt='(F7.5,F15.6,F17.12)') x(i), y(j), Q(1,i,j)
+        end do
+     end do
+     write(1,*)
+     
+     write(1,fmt='(A,I12,I12,I12)') 'VERTICES',n,n*(m+1)
+     do j=1,n
+        write(1,fmt='(I12)', advance='no') m
+        do i=0,m-1
+           write(1,fmt='(I12)', advance='no') i+(j-1)*(m)
+        end do
+        write(1,*)
+     end do
+
+     !
+     ! Store lookup table
+     !
+     write(1,fmt='(A,I12)') 'POINT_DATA',m*n
+     write(1,fmt='(A)') 'SCALARS height double 1'
+     write(1,fmt='(A)') 'LOOKUP_TABLE default'
+     do j=1,n
+        do i=1,m
+           write(1,fmt='(F15.12)') Q(1,i,j)
+        end do
+     end do
+     write(1,*)
+     close(1)
+     
+   end subroutine save_vtk
+end module vtk_export

advanced_lab/run.csh

+#!/bin/csh
+foreach n (`seq 1 1 16`)
+    env OMP_NUM_THREADS=$n aprun -n 1 -d $n ./shwater2d
+end