Merge pull request #626 from nmizukami/advection-duffsion_improved

advection-diffusion equation solver

Author: nmizukami

route/build/src/advection_diffusion.f90

@@ -1,6 +1,6 @@
 MODULE advection_diffusion
 
-! solving advection-diffusion equation
+! Numerical solution for 1D advection-diffusion equation
 
 USE nrtype
 USE public_var,    ONLY: iulog           ! i/o logical unit number
@@ -14,7 +14,7 @@ MODULE advection_diffusion
 CONTAINS
 
   ! *********************************************************************
-  ! subroutine: solve advection-diffuision equation (ade)
+  ! subroutine: solve advection-diffuision equation (ade) per a reach
   ! *********************************************************************
   SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
                        nMolecule,     & ! input: number of sub-segments
@@ -26,23 +26,30 @@ SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
                        FluxPrev,      & ! input: Flux at previous time step [unit of quantity]
                        FluxLocal,     & ! inout: Flux soloved at current time step [unit of quantity]
                        verbose,       & ! input: reach index to be examined
-                       ierr,message)
+                       advec_scheme,  & ! optional input: advection term descretization: 1->central difference, 2->upwind method
+                       downstreamBC,  & ! optional input: downstream end B.C. 1->Neumann, 2->absorbing
+                       wc,            & ! optional input: advection term descretization weight: 0=fully explict to 1=fully implicit
+                       wd             & ! optional input: diffusion term descretization weight: 1=fully explict to 2=fully implicit
+                       )
   ! ----------------------------------------------------------------------------------------
   ! Solve linearlized advection diffusive equation per reach and time step.
   !  dQ/dt + ck*dQ/dx = dk*d2Q/dx2 + FluxLat - a)
   !
   !  ck (velocity) and dk (diffusivity) are given by input argument
   !
   !  1) dQ/dt   = (Q(t,x) - Q(t-1,x))/dt
-  !  2) dQ/dx   = [(1-wck)(Q(t-1,x+1)-Q(t-1,x-1)) + wck*(Q(t,x+1)-Q(t,x-1))]/2dx  (using central difference)
+  !  2) dQ/dx   = [(1-wck)(Q(t-1,x+1)-Q(t-1,x-1)) + wck*(Q(t,x+1)-Q(t,x-1))]/2dx  (central difference)
+  !               or
+  !               [(1-wck)(Q(t-1,x+1)-Q(t-1,x-1)) + wck*(Q(t,x+1)-Q(t,x-1))]/2dx  (upwind)
   !  3) d2Q/dx2 = [(1-wdk)(Q(t-1,x+1)-2Q(t-1,x)+Q(t-1,x-1)) + wdk*(Q(t,x+1)-2Q(t,x)+Q(t,x-1))]/2dx
   !
-  !  upstream B.C:   Dirchlet BC with inflow at current time-step,t, from upstream basin
-  !  downstream B.C: Neumann BC with prescribed quantity gradient (Sbc)
+  !  4) downstream B.C: Neumann BC with prescribed quantity gradient (Sbc)
   !                  dQ/dx|x=N = Sbc ->  4) Q(t,N)-Q(t,N-1)) = Sbc*dx
-  !  Another downstream B.C option is absorbing boundary condition
+  !     Another downstream B.C option: absorbing boundary condition
   !                  dQ/dt|x=N + ck*dQ/dx|x=N = 0
   !
+  !  upstream B.C:   Dirchlet BC with inflow at current time-step,t, from upstream basin
+  !
   !  Inserting 1), 2), 3) and 4) into a) and moving Q(t,*) terms to left-hand side and Q(t-1,*) terms to the righ-hand side
   !  results in tridiagonal matrix equation A*Q = b
   !  where A is [N x N] matrix, Q is [N x 1] vector to be solved (next time step Q) and b is [N x 1] vector
@@ -65,8 +72,10 @@ SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
   real(dp),     intent(in)      :: FluxPrev(nMolecule)       ! sub-reach quantity at previous time step [m3/s]
   real(dp),     intent(out)     :: FluxLocal(nMolecule,0:1)  ! sub-reach & sub-time step quantity at previous and current time step [m3/s]
   logical(lgt), intent(in)      :: verbose                   ! reach index to be examined
-  integer(i4b), intent(out)     :: ierr                      ! error code
-  character(*), intent(out)     :: message                   ! error message
+  integer(i4b), optional, intent(in) :: advec_scheme         ! advection term descretization: 1->central diff.(default), 2->upwind method
+  integer(i4b), optional, intent(in) :: downstreamBC         ! downstream end B.C. 1->Neumann (default), 2->absorbing
+  real(dp), optional, intent(in) :: wc                       ! advection term weight
+  real(dp), optional, intent(in) :: wd                       ! diffusion term weight
   ! Local variables
   real(dp)                      :: Cd                        ! Fourier number
   real(dp)                      :: Ca                        ! Courant number
@@ -79,18 +88,31 @@ SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
   real(dp)                      :: wdk                       ! weight for diffusion
   integer(i4b)                  :: ix                        ! loop index
   integer(i4b)                  :: Nx                        ! number of internal reach segments
-  integer(i4b)                  :: downstreamBC              ! method of B.C condition - absorbing or Neumann
+  integer(i4b)                  :: downBC                    ! B.C condition method used
+  integer(i4b)                  :: advec_discretization      ! advection term discretization used
   character(len=strLen)         :: fmt1                      ! format string
   ! Local parameters
   integer(i4b), parameter       :: absorbingBC=1             ! downstream B.C.
   integer(i4b), parameter       :: neumannBC=2               ! downstream B.C. flux derivative w.r.t. distance at downstream boundary
+  integer(i4b), parameter       :: upwind=1                  ! upwind method for advection term discretization
+  integer(i4b), parameter       :: central=2                 ! central difference for advetion term discretization
 
-  ierr=0; message='solve_ade/'
-
-  ! hard-coded parameters
-  downstreamBC = neumannBC  ! downstream boundary condition
-  wck = 1.0                 ! weight in advection term
-  wdk = 1.0                 ! weight in diffusion term 0.0-> fully explicit, 0.5-> Crank-Nicolson, 1.0 -> fully implicit
+  downBC = neumannBC  ! downstream boundary condition
+  if (present(downstreamBC)) then
+    downBC = downstreamBC
+  end if
+  advec_discretization = central ! advection term discretization method
+  if (present(advec_scheme)) then
+    advec_discretization = advec_scheme
+  end if
+  wck = 1.0   ! weight in advection term: 0.0->fully explict, 1.0->fully implicit (default)
+  if (present(wc)) then
+    wck=wc
+  end if
+  wdk = 1.0   ! weight in diffusion term: 0.0->fully explicit, 0.5->Crank-Nicolson, 1.0->fully implicit (default)
+  if (present(wd)) then
+    wdk = wd
+  end if
 
   Nx = nMolecule - 1  ! Nx: number of internal reach segments
 
@@ -107,47 +129,66 @@ SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
   FluxLocal(1,1)  = FluxUpstream          ! quantity from upstream at current time step
 
   ! Fourier number and Courant number
-  Cd = dk*dt_local/dx**2
+  Cd = dk*dt_local/(dx*dx)
   Ca = ck*dt_local/dx
 
   ! create a matrix - current time step
   ! populate tridiagonal elements
   ! diagonal
   diagonal(1,2)             = 1._dp
-  diagonal(2:nMolecule-1,2) = 2._dp + 4*wdk*Cd
-  if (downstreamBC == absorbingBC) then
+  if (advec_discretization==upwind) then
+    diagonal(2:nMolecule-1,2) = 1._dp + wck*Ca + 2._dp*wdk*Cd
+  else if (advec_discretization==central) then
+    diagonal(2:nMolecule-1,2) = 2._dp + 4*wdk*Cd
+  end if
+  if (downBC == absorbingBC) then
     diagonal(nMolecule,2)     = 1._dp + wck*Ca
-  else if (downstreamBC == neumannBC) then
+  else if (downBC == neumannBC) then
     diagonal(nMolecule,2)     = 1._dp
   end if
 
   ! upper
   diagonal(:,1)           = 0._dp
-  diagonal(3:nMolecule,1) = wck*Ca - 2._dp*wdk*Cd
+  if (advec_discretization==upwind) then
+    diagonal(3:nMolecule,1) = -wdk*Cd
+  else if (advec_discretization==central) then
+    diagonal(3:nMolecule,1) = wck*Ca - 2._dp*wdk*Cd
+  end if
 
   ! lower
   diagonal(:,3)             = 0._dp
-  diagonal(1:nMolecule-2,3) = -wck*Ca - 2._dp*wdk*Cd
-  if (downstreamBC == absorbingBC) then
+  if (advec_discretization==upwind) then
+    diagonal(1:nMolecule-2,3) = -wck*Ca - wdk*Cd
+  else if (advec_discretization==central) then
+    diagonal(1:nMolecule-2,3) = -wck*Ca - 2._dp*wdk*Cd
+  end if
+  if (downBC == absorbingBC) then
     diagonal(nMolecule-1,3)   = -wck*Ca
-  else if (downstreamBC == neumannBC) then
+  else if (downBC == neumannBC) then
     diagonal(nMolecule-1,3)     = -1._dp
   end if
 
   ! populate right-hand side
   ! upstream boundary condition
   b(1)             = FluxLocal(1,1)
   ! downstream boundary condition
-  if (downstreamBC == absorbingBC) then
+  if (downBC == absorbingBC) then
     b(nMolecule) = (1._dp-(1._dp-wck)*Ca)*FluxLocal(nMolecule,0) + (1-wck)*Ca*FluxLocal(nMolecule-1,0)
-  else if (downstreamBC == neumannBC) then
+  else if (downBC == neumannBC) then
     Sbc = (FluxLocal(nMolecule,0)-FluxLocal(nMolecule-1,0))
     b(nMolecule)     = Sbc
   end if
   ! internal node points
-  b(2:nMolecule-1) = ((1._dp-wck)*Ca +2._dp*(1._dp-wdk)*Cd)*FluxLocal(1:nMolecule-2,0)  &
-                    + (2._dp-4._dp*(1._dp-wdk)*Cd)*FluxLocal(2:nMolecule-1,0)           &
-                    - ((1._dp-wck)*Ca -2._dp*(1._dp-wdk)*Cd)*FluxLocal(3:nMolecule,0)
+  if (advec_discretization==upwind) then
+    b(2:nMolecule-1) = ((1._dp-wck)*Ca +(1._dp-wdk)*Cd)*FluxLocal(1:nMolecule-2,0)  &
+                     + (1._dp-(1._dp-wck)*Ca-2._dp*(1._dp-wdk)*Cd)*FluxLocal(2:nMolecule-1,0) &
+                     + (1._dp-wdk)*Cd*FluxLocal(3:nMolecule,0)
+  else if (advec_discretization==central) then
+    b(2:nMolecule-1) = ((1._dp-wck)*Ca +2._dp*(1._dp-wdk)*Cd)*FluxLocal(1:nMolecule-2,0)  &
+                      + (2._dp-4._dp*(1._dp-wdk)*Cd)*FluxLocal(2:nMolecule-1,0)           &
+                      - ((1._dp-wck)*Ca -2._dp*(1._dp-wdk)*Cd)*FluxLocal(3:nMolecule,0)
+  end if
+
   ! solve matrix equation - get updated FluxLocal
   call TDMA(nMolecule, diagonal, b, FluxSolved)
 
@@ -162,7 +203,7 @@ SUBROUTINE solve_ade(reach_length,  & ! input: Reach length [m]
     write(iulog,fmt1) ' diagonal(:,2)= ', diagonal(:,2)
     write(iulog,fmt1) ' diagonal(:,3)= ', diagonal(:,3)
     write(iulog,fmt1) ' b= ', b(1:nMolecule)
-    write(iulog,fmt1) ' Q sub_reqch=', (FluxSolved(ix), ix=1,nMolecule)
+    write(iulog,fmt1) ' FluxSolved= ', (FluxSolved(ix), ix=1,nMolecule)
   end if
 
   END SUBROUTINE solve_ade

route/build/src/dfw_route.f90

@@ -33,6 +33,8 @@ MODULE dfw_route_module
 
 integer(i4b), parameter :: top_reach=1
 integer(i4b), parameter :: bottom_reach=2
+integer(i4b), parameter :: advec_scheme=2 ! 1->upwind, 2->central difference
+integer(i4b), parameter :: downBC = 2 ! downstream end B.C: 1->open (absorbing) B.C., 2->Neumann B.C.
 
 type, extends(base_route_rch) :: dfw_route_rch
  CONTAINS
@@ -253,7 +255,8 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
 
  ntSub = 1  ! number of sub-time step
 
- associate(S         => rch_param%R_SLOPE,    & ! channel slope
+ associate(nMolecule => nMolecule%DW_ROUTE,   & ! number of computing points in a reach
+           S         => rch_param%R_SLOPE,    & ! channel slope
            n         => rch_param%R_MAN_N,    & ! manning n
            bt        => rch_param%R_WIDTH,    & ! channel bottom width
            bankDepth => rch_param%R_DEPTH,    & ! bankfull depth
@@ -266,11 +269,11 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
 
    if (L > min_length_route) then
 
-   allocate(Qprev(nMolecule%DW_ROUTE), stat=ierr, errmsg=cmessage)
+   allocate(Qprev(nMolecule), stat=ierr, errmsg=cmessage)
    if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
 
    ! initialize previous time step flow
-   Qprev(1:nMolecule%DW_ROUTE) = rstate%molecule%Q     ! flow state at previous time step
+   Qprev(1:nMolecule) = rstate%molecule%Q     ! flow state at previous time step
 
    if (verbose) then
      write(iulog,'(A,1X,G12.5)') ' length [m]        =',L
@@ -286,17 +289,17 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
      write(iulog,'(A,1X,I3,A,1X,G12.5)') ' No. sub timestep=',nTsub,' sub time-step [sec]=',dTsub
    end if
 
-   allocate(Qlocal(1:nMolecule%DW_ROUTE, 0:1), stat=ierr, errmsg=cmessage)
+   allocate(Qlocal(1:nMolecule, 0:1), stat=ierr, errmsg=cmessage)
    if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
 
    do it = 1, nTsub
-     Qbar = (Qupstream+Qprev(1)+Qprev(nMolecule%DW_ROUTE-1))/3.0 ! 3 point average discharge [m3/s]
+     Qbar = (Qupstream+Qprev(1)+Qprev(nMolecule-1))/3.0 ! 3 point average discharge [m3/s]
      depth = flow_depth(abs(Qbar), bt, zc, S, n, zf=zf, bankDepth=bankDepth) ! compute flow depth as normal depth (a function of flow)
      ck    = celerity(abs(Qbar), depth, bt, zc, S, n, zf=zf, bankDepth=bankDepth)
      dk    = diffusivity(abs(Qbar), depth, bt, zc, S, n, zf=zf, bankDepth=bankDepth)
 
      call solve_ade(L,                  & ! input: river parameter data structure
-                    nMolecule%DW_ROUTE, & ! input: number of sub-segments
+                    nMolecule,          & ! input: number of sub-segments
                     dTsub,              & ! input: time_step [sec]
                     Qupstream,          & ! input: quantity from upstream [unit of quantity]
                     ck,                 & ! input: velocity [m/s]
@@ -305,19 +308,19 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
                     Qprev,              & ! input: quantity at previous time step [unit of quantity]
                     Qlocal,             & ! inout: quantity soloved at current time step [unit of quantity]
                     verbose,            & ! input: reach index to be examined
-                    ierr,message)
-     if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
+                    advec_scheme=advec_scheme, & ! optional input: advection scheme (1->central difference, 2->upwind)
+                    downstreamBC=downBC)         ! optional input: downstream boundary condition (1->Neumann, 2->open boundary condition)
    end do
 
    ! For very low flow condition, outflow - inflow may exceed current storage, so limit outflow and adjust flow profile
-   if (abs(Qlocal(nMolecule%DW_ROUTE-1,1))>0._dp) then
+   if (abs(Qlocal(nMolecule-1,1))>0._dp) then
      volTmp = max(0._dp, rflux%ROUTE(idxDW)%REACH_VOL(1))
-     qoutTmp = Qlocal(nMolecule%DW_ROUTE-1,1)*dt
+     qoutTmp = Qlocal(nMolecule-1,1)*dt
      pcntReduc = min((volTmp + dt*Qupstream)*0.999_dp/qoutTmp, 1._dp)
-     Qlocal(2:nMolecule%DW_ROUTE,1) = Qlocal(2:nMolecule%DW_ROUTE,1)*pcntReduc
+     Qlocal(2:nMolecule,1) = Qlocal(2:nMolecule,1)*pcntReduc
    end if
 
-   rflux%ROUTE(idxDW)%REACH_VOL(1) = rflux%ROUTE(idxDW)%REACH_VOL(1) + (Qupstream - Qlocal(nMolecule%DW_ROUTE-1,1))*dt
+   rflux%ROUTE(idxDW)%REACH_VOL(1) = rflux%ROUTE(idxDW)%REACH_VOL(1) + (Qupstream - Qlocal(nMolecule-1,1))*dt
 
    ! if reach volume exceeds flood threshold volume, excess water is flooded volume.
    if (rflux%ROUTE(idxDW)%REACH_VOL(1) > bankVol) then
@@ -329,15 +332,15 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
    rflux%ROUTE(idxDW)%REACH_ELE = water_height(rflux%ROUTE(idxDW)%REACH_VOL(1)/L, bt, zc, zf=zf, bankDepth=bankDepth)
 
    ! store final outflow in data structure
-   rflux%ROUTE(idxDW)%REACH_Q = Qlocal(nMolecule%DW_ROUTE-1,1) + Qlat
+   rflux%ROUTE(idxDW)%REACH_Q = Qlocal(nMolecule-1,1) + Qlat
 
    ! update state
    rstate%molecule%Q = Qlocal(:,1)
 
    else ! length < min_length_route: length is short enough to just pass upstream to downstream
      rflux%ROUTE(idxDW)%REACH_Q = Qupstream + Qlat
-     rstate%molecule%Q(1:nMolecule%DW_ROUTE) = 0._dp
-     rstate%molecule%Q(nMolecule%DW_ROUTE)   = rflux%ROUTE(idxDW)%REACH_Q
+     rstate%molecule%Q(1:nMolecule) = 0._dp
+     rstate%molecule%Q(nMolecule)   = rflux%ROUTE(idxDW)%REACH_Q
 
      rflux%ROUTE(idxDW)%REACH_VOL(0) = 0._dp
      rflux%ROUTE(idxDW)%REACH_VOL(1) = 0._dp
@@ -353,8 +356,8 @@ SUBROUTINE diffusive_wave(rch_param,     & ! input: river parameter data structu
    rflux%ROUTE(idxDW)%FLOOD_VOL(1) = 0._dp
    rflux%ROUTE(idxDW)%REACH_ELE    = 0._dp
 
-   rstate%molecule%Q(1:nMolecule%DW_ROUTE) = 0._dp
-   rstate%molecule%Q(nMolecule%DW_ROUTE)   = rflux%ROUTE(idxDW)%REACH_Q
+   rstate%molecule%Q(1:nMolecule) = 0._dp
+   rstate%molecule%Q(nMolecule)   = rflux%ROUTE(idxDW)%REACH_Q
 
    if (verbose) then
      write(iulog,'(A)')            ' This is headwater '

route/build/src/kwe_route.f90

@@ -302,9 +302,7 @@ SUBROUTINE kinematic_wave(rch_param,     & ! input: river parameter data structu
                     Qlat,               & ! input: lateral quantity into chaneel [unit of quantity]
                     Qprev,              & ! input: quantity at previous time step [unit of quantity]
                     Qlocal,             & ! inout: quantity soloved at current time step [unit of quantity]
-                    verbose,            & ! input: reach index to be examined
-                    ierr,message)
-     if(ierr/=0)then; message=trim(message)//trim(cmessage); return; endif
+                    verbose)              ! input: reach index to be examined
    end do
 
    ! For very low flow condition, outflow - inflow may exceed current storage, so limit outflow and adjust flow profile