Add linear drag and use impulse on wall for kelvin wave setup

Author: fluidnumericsJoe

examples/linear_shallow_water2d_kelvinwaves.f90

@@ -33,9 +33,10 @@ program linear_shallow_water2d_kelvinwaves
   character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' ! Which integrator method
   integer,parameter :: controlDegree = 7 ! Degree of control polynomial
   integer,parameter :: targetDegree = 16 ! Degree of target polynomial
-  real(prec),parameter :: dt = 0.001_prec ! Time-step size
-  real(prec),parameter :: endtime = 1.0_prec !30.0_prec ! (s);
-  real(prec),parameter :: f0 = -10.0_prec ! reference coriolis parameter (1/s)
+  real(prec),parameter :: dt = 0.0025_prec ! Time-step size
+  real(prec),parameter :: endtime = 30.0_prec ! (s);
+  real(prec),parameter :: f0 = 10.0_prec ! reference coriolis parameter (1/s)
+  real(prec),parameter :: Cd = 0.5_prec ! Linear drag coefficient (1/s)
   real(prec),parameter :: iointerval = 0.05 ! Write files 20 times per characteristic time scale
   real(prec) :: r
   real(prec) :: e0,ef ! Initial and final entropy
@@ -71,25 +72,26 @@ program linear_shallow_water2d_kelvinwaves
   ! Set the resting surface height and gravity
   modelobj%H = H
   modelobj%g = g
+  modelobj%Cd = Cd
 
   ! ! Set the initial conditions
-  !call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( x^2 + y^2 )/0.02 ) ')
-  !call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
+  call modelobj%solution%SetEquation(3,'f = 0.001*exp( -( (x-1.0)^2 + y^2 )/0.02 ) ')
+  call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec)
 
-  do iel = 1,modelobj%mesh%nElem
-    do j = 1,modelobj%solution%N+1
-      do i = 1,modelobj%solution%N+1
-        call random_number(r)
-        !  modelobj%solution%interior(i,j,iel,3) = modelobj%solution%interior(i,j,iel,3) + 0.0001_prec*(r-0.5)
-        modelobj%solution%interior(i,j,iel,3) = 0.0001_prec*(r-0.5)
+  ! do iel = 1,modelobj%mesh%nElem
+  !   do j = 1,modelobj%solution%N+1
+  !     do i = 1,modelobj%solution%N+1
+  !       call random_number(r)
+  !       modelobj%solution%interior(i,j,iel,3) = modelobj%solution%interior(i,j,iel,3)*(1.0_prec+1.0_prec*(r-0.5))
+  !      ! modelobj%solution%interior(i,j,iel,3) = 0.0001_prec*(r-0.5)
 
-      enddo
-    enddo
-  enddo
-  call modelobj%solution%UpdateDevice()
+  !     enddo
+  !   enddo
+  ! enddo
+  ! call modelobj%solution%UpdateDevice()
 
   call modelobj%SetCoriolis(f0)
-  call modelobj%DiagnoseGeostrophicVelocity()
+  ! call modelobj%DiagnoseGeostrophicVelocity()
 
   call modelobj%WriteModel()
   call modelobj%IncrementIOCounter()
@@ -108,9 +110,8 @@ program linear_shallow_water2d_kelvinwaves
 
   ef = modelobj%entropy
 
-  print*,e0,ef
-  if(abs(ef-e0) > epsilon(e0)) then
-    print*,"Warning: Final entropy greater than initial entropy! ",e0,ef
+  if(ef > e0) then
+    print*,"Warning: Final entropy greater than initial entropy! ",ef,e0
   endif
 
   ! Clean up

src/SELF_LinearShallowWater2D_t.f90

@@ -34,6 +34,7 @@ module self_LinearShallowWater2D_t
   type,extends(dgmodel2d) :: LinearShallowWater2D_t
     real(prec) :: H = 0.0_prec ! uniform resting depth
     real(prec) :: g = 0.0_prec ! acceleration due to gravity
+    real(prec) :: Cd = 0.0_prec ! Linear drag coefficient (1/s)
     type(MappedScalar2D) :: fCori ! The coriolis parameter
 
   contains
@@ -246,8 +247,8 @@ subroutine sourcemethod_LinearShallowWater2D_t(this)
 
       s = this%solution%interior(i,j,iel,1:this%nvar)
 
-      this%source%interior(i,j,iel,1) = this%fCori%interior(i,j,iel,1)*s(2) ! du/dt = f*v
-      this%source%interior(i,j,iel,2) = -this%fCori%interior(i,j,iel,1)*s(1) ! dv/dt = -f*u
+      this%source%interior(i,j,iel,1) = this%fCori%interior(i,j,iel,1)*s(2)-this%Cd*s(1) ! du/dt = f*v - Cd*u
+      this%source%interior(i,j,iel,2) = -this%fCori%interior(i,j,iel,1)*s(1)-this%Cd*s(2) ! dv/dt = -f*u - Cd*v
 
     enddo
 

src/gpu/SELF_LinearShallowWater2D.cpp

@@ -144,25 +144,25 @@ extern "C"
   }
 }
 
-__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, int ndof){
+__global__ void sourcemethod_LinearShallowWater2D_gpukernel(real *solution, real *source, real *fCori, real Cd, int ndof){
   uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x;
 
   if( idof < ndof ){
     real u = solution[idof];
     real v = solution[idof + ndof];
 
-    source[idof] = fCori[idof]*v; // du/dt = fv
-    source[idof+ndof] = -fCori[idof]*u;   // dv/dt  = -fu
+    source[idof] = fCori[idof]*v - Cd*u; // du/dt = fv - Cd*u
+    source[idof+ndof] = -fCori[idof]*u - Cd*v;   // dv/dt  = -fu - Cd*v
 
   }
 
 }
 extern "C"
 {
-  void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, int N, int nel, int nvar){
+  void sourcemethod_LinearShallowWater2D_gpu(real *solution, real *source, real *fCori, real Cd, int N, int nel, int nvar){
     int ndof = (N+1)*(N+1)*nel;
     int threads_per_block = 256;
     int nblocks_x = ndof/threads_per_block +1;
-    sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,ndof);
+    sourcemethod_LinearShallowWater2D_gpukernel<<<dim3(nblocks_x,1,1), dim3(threads_per_block,1,1), 0, 0>>>(solution,source,fCori,Cd,ndof);
   }
 }

src/gpu/SELF_LinearShallowWater2D.f90

@@ -71,11 +71,12 @@ subroutine fluxmethod_LinearShallowWater2D_gpu(solution,flux,g,H,N,nel,nvar) &
   endinterface
 
   interface
-    subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,N,nel,nvar) &
+    subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,Cd,N,nel,nvar) &
       bind(c,name="sourcemethod_LinearShallowWater2D_gpu")
       use iso_c_binding
       use SELF_Constants
       type(c_ptr),value :: solution,source,fCori
+      real(c_prec),value :: Cd
       integer(c_int),value :: N,nel,nvar
     endsubroutine sourcemethod_LinearShallowWater2D_gpu
   endinterface
@@ -165,6 +166,7 @@ subroutine sourcemethod_LinearShallowWater2D(this)
     call sourcemethod_LinearShallowWater2D_gpu(this%solution%interior_gpu, &
                                                this%source%interior_gpu, &
                                                this%fCori%interior_gpu, &
+                                               this%Cd, &
                                                this%solution%interp%N, &
                                                this%solution%nelem, &
                                                this%solution%nvar)