diff --git a/src/__init__.py b/src/__init__.py index fbcdbbed..739c1950 100755 --- a/src/__init__.py +++ b/src/__init__.py @@ -63,6 +63,7 @@ from . import nonlinear_elasticity_fwave_1D from . import reactive_euler_with_efix_1D from . import shallow_roe_with_efix_1D + from . import shallow_bathymetry_fwave_1D from . import shallow_roe_tracer_1D from . import traffic_1D from . import traffic_vc_1D @@ -79,6 +80,8 @@ from . import kpp_2D from . import psystem_2D from . import shallow_roe_with_efix_2D + from . import shallow_bathymetry_fwave_2D + from . import sw_aug_2D from . import shallow_sphere_2D from . import vc_acoustics_2D from . import vc_advection_2D diff --git a/src/rp1_shallow_bathymetry_fwave.f90 b/src/rp1_shallow_bathymetry_fwave.f90 new file mode 100644 index 00000000..3a8355b8 --- /dev/null +++ b/src/rp1_shallow_bathymetry_fwave.f90 @@ -0,0 +1,113 @@ +subroutine rp1(maxm,num_eqn,num_waves,num_aux,num_ghost,num_cells,ql,qr,auxl,auxr,fwave,s,amdq,apdq) + +! Riemann solver for the 1D shallow water equations: + +! (h)_t + (h*u)_x = 0 +! (hu)_t + (h*u^2 + 1/2*grav*h^2)_x = -grav*h*(b)_x +! +! using the f-wave algorithm and Roe's approximate Riemann solver. + +! waves: 2 +! equations: 2 + +! Conserved quantities: +! 1 depth +! 2 momentum + +! Auxiliary fields: +! 1 bathymetry + +! The gravitational constant grav should be in the common block cparam. +! +! See http://www.clawpack.org/riemann.html for a detailed explanation +! of the Riemann solver API. + + implicit none + + integer, intent(in) :: maxm, num_eqn, num_waves, num_ghost, num_aux, num_cells + real(kind=8), intent(in) :: ql(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: qr(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxl(num_aux, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxr(num_aux, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: s(num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: fwave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: amdq(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: apdq(num_eqn, 1-num_ghost:maxm+num_ghost) + + real(kind=8) :: grav, dry_tolerance, sea_level + common /cparam/ grav, dry_tolerance, sea_level + + real(kind=8) :: hl, ul, hr, ur, hbar, uhat, chat, bl, br, phil, phir + real(kind=8) :: R(2,2), fluxdiff(2), beta(2) + + integer :: i, k + + amdq = 0.d0 + apdq = 0.d0 + + do i=2-num_ghost,num_cells+num_ghost + ! # Left states + hl = qr(1, i - 1) + if (hl > dry_tolerance) then + ul = qr(2, i - 1) / hl + phil = qr(1, i - 1) * ul**2 + 0.5d0 * grav * qr(1, i - 1)**2 + else + ul = 0.d0 + phil = 0.d0 + end if + bl = auxr(1, i - 1) + + ! # Right states + hr = ql(1, i) + if (hr > dry_tolerance) then + ur = ql(2, i)/hr + phir = ql(1, i) * ur**2 + 0.5d0 * grav * ql(1, i)**2 + else + ur = 0.d0 + phir = 0.d0 + end if + br = auxl(1, i) + + ! # Roe average states + hbar = 0.5 * (hr + hl) + uhat = (sqrt(hl) * ul + sqrt(hr) * ur) / (sqrt(hl) + sqrt(hr)) + chat = sqrt(grav * hbar) + + ! # Flux differences + fluxdiff(1) = (hr * ur) - (hl * ul) + fluxdiff(2) = phir - phil + grav * hbar * (br - bl) + + ! # Wave speeds + s(1,i) = min(ul - sqrt(grav * hl), uhat - chat) + s(2,i) = max(ur + sqrt(grav * hr), uhat + chat) + + ! # Right eigenvectors (column) + R(1,1) = 1.d0 + R(2,1) = s(1, i) + + R(1,2) = 1.d0 + R(2,2) = s(2, i) + + ! Wave strengths + beta(1) = (s(2, i) * fluxdiff(1) - fluxdiff(2)) / (s(2, i) - s(1, i)) + beta(2) = (fluxdiff(2) - s(1, i) * fluxdiff(1)) / (s(2, i) - s(1, i)) + + ! # Flux waves + do k=1,num_waves + fwave(:,k,i) = beta(k) * R(:,k) + enddo + + ! # Fluctuations + do k=1, num_waves + if (s(k, i) < 0.d0) then + amdq(:, i) = amdq(:, i) + fwave(:, k, i) + else if (s(k, i) > 0.d0) then + apdq(:, i) = apdq(:, i) + fwave(:, k, i) + else + amdq(:, i) = amdq(:, i) + 0.5d0 * fwave(:, k, i) + apdq(:, i) = apdq(:, i) + 0.5d0 * fwave(:, k, i) + end if + enddo + enddo + +end subroutine rp1 diff --git a/src/rpn2_shallow_bathymetry_fwave.f90 b/src/rpn2_shallow_bathymetry_fwave.f90 new file mode 100644 index 00000000..a59555ba --- /dev/null +++ b/src/rpn2_shallow_bathymetry_fwave.f90 @@ -0,0 +1,142 @@ +subroutine rpn2(ixy, maxm, num_eqn, num_waves, num_aux, num_ghost, & + num_cells, ql, qr, auxl, auxr, fwave, s, amdq, apdq) + +! Riemann solver for the 2D shallow water equations + +! (h)_t + (h*u)_x = 0 +! (hu)_t + (h*u^2 + 1/2*grav*h^2)_x + (h*u*v)_y = -grav*h*(b)_x +! (hv)_t + (h*u*v)_x + (h*v^2 + 1/2*grav*h^2)_y = -grav*h*(b)_y + +! using the f-wave algorithm and Roe's approximate Riemann solver. + +! waves: 3 +! equations: 3 + +! Conserved quantities: +! 1 depth +! 2 x_momentum +! 3 y_momentum + +! Auxiliary fields: +! 1 bathymetry + +! The gravitational constant grav should be in the common block cparam. +! +! See http://www.clawpack.org/riemann.html for a detailed explanation +! of the Riemann solver API. + + implicit none + + integer, intent(in) :: ixy, maxm, num_eqn, num_waves, num_ghost, num_aux, num_cells + real(kind=8), intent(in) :: ql(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: qr(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxl(num_aux, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxr(num_aux, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: s(num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: fwave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: amdq(num_eqn,1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: apdq(num_eqn,1-num_ghost:maxm+num_ghost) + + real(kind=8) :: grav, dry_tolerance, sea_level + common /cparam/ grav, dry_tolerance, sea_level + + real(kind=8) :: hl, ul, vl, bl, hr, ur, vr, br, hbar, uhat, chat + real(kind=8) :: phil, phir, dry_state_l, dry_state_r + real(kind=8) :: R(3,3) + real(kind=8) :: fluxdiff(3), beta(3) + + integer :: i, k, normal_index, transverse_index + + ! Determine normal and tangential directions + if (ixy == 1) then + normal_index = 2 + transverse_index = 3 + else + normal_index = 3 + transverse_index = 2 + end if + + amdq = 0.d0 + apdq = 0.d0 + + ! Primary loop over each cell + do i = 2 - num_ghost, num_cells + num_ghost + + ! Check for dry states - need merge here to convert to float + dry_state_l = merge(0.d0, 1.d0, qr(1, i - 1) < dry_tolerance) + dry_state_r = merge(0.d0, 1.d0, ql(1, i) < dry_tolerance) + + ! Note that for the states below u is always the normal velocity and + ! v is always the tangential velocity + + ! Left states + hl = qr(1, i - 1) * dry_state_l + ul = qr(normal_index, i - 1) / qr(1, i - 1) * dry_state_l + vl = qr(transverse_index, i - 1) / qr(1, i - 1) * dry_state_l + phil = (0.5d0 * grav * hl**2 + hl * ul**2) * dry_state_l + + bl = auxr(1, i - 1) + + ! Right states + hr = ql(1, i) * dry_state_r + ur = ql(normal_index, i) / ql(1, i) * dry_state_r + vr = ql(transverse_index, i) / ql(1, i) * dry_state_r + phir = (0.5d0 * grav * hr**2 + hr * ur**2) * dry_state_r + + br = auxl(1, i) + + ! Roe average states (Roe's linearization) + hbar = 0.5d0 * (hr + hl) + uhat = (sqrt(hr) * ur + sqrt(hl) * ul) / (sqrt(hr) + sqrt(hl)) + chat = sqrt(grav * hbar) + + ! Flux differences + fluxdiff(1) = hr * ur - hl * ul + fluxdiff(2) = phir - phil + grav * hbar * (br - bl) + fluxdiff(3) = hr * ur * vr - hl * ul * vl + + ! Wave speeds + s(1, i) = min(uhat - chat, ul - sqrt(grav * hl)) + s(3, i) = max(uhat + chat, ur + sqrt(grav * hr)) + s(2, i) = 0.5d0 * (s(1, i) + s(3, i)) + + ! Right eigenvectors (columns) + ! could possibly use vhat instead of vl and vr + R(1, 1) = 1.d0 + R(normal_index, 1) = s(1, i) + R(transverse_index, 1) = vl + + R(1, 2) = 0.d0 + R(normal_index, 2) = 0.0 + R(transverse_index, 2) = 1.0 + + R(1, 3) = 1.d0 + R(normal_index, 3) = s(3, i) + R(transverse_index, 3) = vr + + ! Wave strengths + beta(1) = (s(3, i) * fluxdiff(1) - fluxdiff(2)) / (s(3, i) - s(1, i)) + beta(3) = (fluxdiff(2) - s(1, i) * fluxdiff(1)) / (s(3, i) - s(1, i)) + beta(2) = fluxdiff(3) - beta(1) * vl - beta(3) * vr + + ! f-waves + do k = 1, num_waves + fwave(:, k, i) = beta(k) * R(:, k) + enddo + + ! Fluctuations - could probably rewrite this to be a masking operation + ! instead... + do k=1, num_waves + if (s(k, i) < 1.0e-14) then + amdq(:, i) = amdq(:, i) + fwave(:, k, i) + elseif (s(k, i) > 1.0e-14) then + apdq(:, i) = apdq(:, i) + fwave(:, k, i) + else + amdq(:, i) = amdq(:, i) + 0.5d0 * fwave(:, k, i) + apdq(:, i) = apdq(:, i) + 0.5d0 * fwave(:, k, i) + endif + enddo + + enddo ! End of main loop + +end subroutine rpn2 diff --git a/src/rpn2_sw_aug.f90 b/src/rpn2_sw_aug.f90 new file mode 100755 index 00000000..9972b218 --- /dev/null +++ b/src/rpn2_sw_aug.f90 @@ -0,0 +1,865 @@ +subroutine rpn2(ixy, maxm, num_eqn, num_waves, num_aux, num_ghost, & + num_cells, ql, qr, auxl, auxr, fwave, s, amdq, apdq) + +! Normal Riemann solver for the 2D SHALLOW WATER equations +! with topography: +! # h_t + (hu)_x + (hv)_y = 0 # +! # (hu)_t + (hu^2 + 0.5gh^2)_x + (huv)_y = -ghb_x # +! # (hv)_t + (huv)_x + (hv^2 + 0.5gh^2)_y = -ghb_y # + +! This solver is based on David George's solver written for GeoClaw. +! It has been modified to be compatible with f2py (and thus PyClaw). + +! waves: 3 +! equations: 3 + +! Conserved quantities: +! 1 depth +! 2 x_momentum +! 3 y_momentum + +! Auxiliary fields: +! 1 bathymetry + +! The gravitational constant grav should be in the common block cparam. + +! See http://www.clawpack.org/riemann.html for a detailed explanation +! of the Riemann solver API. + + implicit none + + real(kind=8) :: grav, g + real(kind=8), parameter :: drytol = 1.e-14 + common /cparam/ grav + + integer, intent(in) :: maxm,num_eqn,num_aux,num_waves,num_ghost,num_cells,ixy + + real(kind=8), intent(inout) :: ql(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(inout) :: qr(num_eqn, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxl(num_aux,1-num_ghost:maxm+num_ghost) + real(kind=8), intent(in) :: auxr(num_aux,1-num_ghost:maxm+num_ghost) + + real(kind=8), intent(out) :: fwave(num_eqn, num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: s(num_waves, 1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: apdq(num_eqn,1-num_ghost:maxm+num_ghost) + real(kind=8), intent(out) :: amdq(num_eqn,1-num_ghost:maxm+num_ghost) + + !local only + integer m,i,mw,maxiter,mu,nv + real(kind=8) wall(3) + real(kind=8) fw(3,3) + real(kind=8) sw(3) + + real(kind=8) hR,hL,huR,huL,uR,uL,hvR,hvL,vR,vL,phiR,phiL + real(kind=8) bR,bL,sL,sR,sRoe1,sRoe2,sE1,sE2,uhat,chat + real(kind=8) s1m,s2m + real(kind=8) hstar,hstartest,hstarHLL,sLtest,sRtest + real(kind=8) tw,dxdc + + logical rare1,rare2 + + g = grav + + !loop through Riemann problems at each grid cell + do i=2-num_ghost,num_cells+num_ghost + + !-----------------------Initializing------------------------------ + !inform of a bad riemann problem from the start + if((qr(1,i-1).lt.0.d0).or.(ql(1,i) .lt. 0.d0)) then + write(*,*) 'Negative input: hl,hr,i=',qr(1,i-1),ql(1,i),i + endif + + !Initialize Riemann problem for grid interface + do mw=1,num_waves + s(mw,i)=0.d0 + fwave(1,mw,i)=0.d0 + fwave(2,mw,i)=0.d0 + fwave(3,mw,i)=0.d0 + enddo + +! !set normal direction + if (ixy.eq.1) then + mu=2 + nv=3 + else + mu=3 + nv=2 + endif + + !zero (small) negative values if they exist + if (qr(1,i-1).lt.0.d0) then + qr(1,i-1)=0.d0 + qr(2,i-1)=0.d0 + qr(3,i-1)=0.d0 + endif + + if (ql(1,i).lt.0.d0) then + ql(1,i)=0.d0 + ql(2,i)=0.d0 + ql(3,i)=0.d0 + endif + + !skip problem if in a completely dry area + if (qr(1,i-1) <= drytol .and. ql(1,i) <= drytol) then + go to 30 + endif + + !Riemann problem variables + hL = qr(1,i-1) + hR = ql(1,i) + huL = qr(mu,i-1) + huR = ql(mu,i) + bL = auxr(1,i-1) + bR = auxl(1,i) + + hvL=qr(nv,i-1) + hvR=ql(nv,i) + + !check for wet/dry boundary + if (hR.gt.drytol) then + uR=huR/hR + vR=hvR/hR + phiR = 0.5d0*g*hR**2 + huR**2/hR + else + hR = 0.d0 + huR = 0.d0 + hvR = 0.d0 + uR = 0.d0 + vR = 0.d0 + phiR = 0.d0 + endif + + if (hL.gt.drytol) then + uL=huL/hL + vL=hvL/hL + phiL = 0.5d0*g*hL**2 + huL**2/hL + else + hL=0.d0 + huL=0.d0 + hvL=0.d0 + uL=0.d0 + vL=0.d0 + phiL = 0.d0 + endif + + wall(1) = 1.d0 + wall(2) = 1.d0 + wall(3) = 1.d0 + if (hR.le.drytol) then + call riemanntype(hL,hL,uL,-uL,hstar,s1m,s2m,rare1,rare2,1,drytol,g) + hstartest=max(hL,hstar) + if (hstartest+bL.lt.bR) then !right state should become ghost values that mirror left for wall problem +! bR=hstartest+bL + wall(2)=0.d0 + wall(3)=0.d0 + hR=hL + huR=-huL + bR=bL + phiR=phiL + uR=-uL + vR=vL + elseif (hL+bL.lt.bR) then + bR=hL+bL + endif + elseif (hL.le.drytol) then ! right surface is lower than left topo + call riemanntype(hR,hR,-uR,uR,hstar,s1m,s2m,rare1,rare2,1,drytol,g) + hstartest=max(hR,hstar) + if (hstartest+bR.lt.bL) then !left state should become ghost values that mirror right +! bL=hstartest+bR + wall(1)=0.d0 + wall(2)=0.d0 + hL=hR + huL=-huR + bL=bR + phiL=phiR + uL=-uR + vL=vR + elseif (hR+bR.lt.bL) then + bL=hR+bR + endif + endif + + !determine wave speeds + sL=uL-sqrt(g*hL) ! 1 wave speed of left state + sR=uR+sqrt(g*hR) ! 2 wave speed of right state + + uhat=(sqrt(g*hL)*uL + sqrt(g*hR)*uR)/(sqrt(g*hR)+sqrt(g*hL)) ! Roe average + chat=sqrt(g*0.5d0*(hR+hL)) ! Roe average + sRoe1=uhat-chat ! Roe wave speed 1 wave + sRoe2=uhat+chat ! Roe wave speed 2 wave + + sE1 = min(sL,sRoe1) ! Eindfeldt speed 1 wave + sE2 = max(sR,sRoe2) ! Eindfeldt speed 2 wave + + !--------------------end initializing...finally---------- + !solve Riemann problem. + + maxiter = 1 + + call riemann_aug_JCP(maxiter,3,3,hL,hR,huL, & + huR,hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2, & + drytol,g,sw,fw) + +! !eliminate ghost fluxes for wall + do mw=1,3 + sw(mw)=sw(mw)*wall(mw) + + fw(1,mw)=fw(1,mw)*wall(mw) + fw(2,mw)=fw(2,mw)*wall(mw) + fw(3,mw)=fw(3,mw)*wall(mw) + enddo + + do mw=1,num_waves + s(mw,i)=sw(mw) + fwave(1,mw,i)=fw(1,mw) + fwave(mu,mw,i)=fw(2,mw) + fwave(nv,mw,i)=fw(3,mw) + enddo + + 30 continue + enddo + + + +!============= compute fluctuations============================================= + amdq(1:3,:) = 0.d0 + apdq(1:3,:) = 0.d0 + do i=2-num_ghost,num_cells+num_ghost + do mw=1,num_waves + if (s(mw,i) < 0.d0) then + amdq(1:3,i) = amdq(1:3,i) + fwave(1:3,mw,i) + else if (s(mw,i) > 0.d0) then + apdq(1:3,i) = apdq(1:3,i) + fwave(1:3,mw,i) + else + amdq(1:3,i) = amdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i) + apdq(1:3,i) = apdq(1:3,i) + 0.5d0 * fwave(1:3,mw,i) + endif + enddo + enddo + + return + end subroutine + + +!----------------------------------------------------------------------- + subroutine riemann_aug_JCP(maxiter,num_eqn,num_waves,hL,hR,huL,huR, & + hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! This solver is described in J. Comput. Phys. (6): 3089-3113, March 2008 + ! Augmented Riemann Solvers for the Shallow Equations with Steady States and Inundation + + ! To use the original solver call with maxiter=1. + + ! This solver allows iteration when maxiter > 1. The iteration seems to help with + ! instabilities that arise (with any solver) as flow becomes transcritical over variable topo + ! due to loss of hyperbolicity. + + + + implicit none + + !input + integer num_eqn,num_waves,maxiter + real(kind=8) fw(num_eqn,num_waves) + real(kind=8) sw(num_waves) + real(kind=8) hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,sE1,sE2 + real(kind=8) hvL,hvR,vL,vR + real(kind=8) drytol,g + + + !local + integer m,mw,k,iter + real(kind=8) A(3,3) + real(kind=8) r(3,3) + real(kind=8) lambda(3) + real(kind=8) del(3) + real(kind=8) beta(3) + + real(kind=8) delh,delhu,delphi,delb,delnorm + real(kind=8) rare1st,rare2st,sdelta,raremin,raremax + real(kind=8) criticaltol,convergencetol,raretol + real(kind=8) s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar + real(kind=8) huRstar,huLstar,uRstar,uLstar,hstarHLL + real(kind=8) deldelh,deldelphi + real(kind=8) s1m,s2m,hm + real(kind=8) det1,det2,det3,determinant + + logical rare1,rare2,rarecorrector,rarecorrectortest,sonic + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + delnorm = delh**2 + delphi**2 + + call riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2,1,drytol,g) + + + lambda(1)= min(sE1,s2m) !Modified Einfeldt speed + lambda(3)= max(sE2,s1m) !Modified Eindfeldt speed + sE1=lambda(1) + sE2=lambda(3) + lambda(2) = 0.d0 ! ### Fix to avoid uninitialized value in loop on mw -- Correct?? ### + + + hstarHLL = max((huL-huR+sE2*hR-sE1*hL)/(sE2-sE1),0.d0) ! middle state in an HLL solve + +! !determine the middle entropy corrector wave------------------------ + rarecorrectortest=.false. + rarecorrector=.false. + if (rarecorrectortest) then + sdelta=lambda(3)-lambda(1) + raremin = 0.5d0 + raremax = 0.9d0 + if (rare1.and.sE1*s1m.lt.0.d0) raremin=0.2d0 + if (rare2.and.sE2*s2m.lt.0.d0) raremin=0.2d0 + if (rare1.or.rare2) then + !see which rarefaction is larger + rare1st=3.d0*(sqrt(g*hL)-sqrt(g*hm)) + rare2st=3.d0*(sqrt(g*hR)-sqrt(g*hm)) + if (max(rare1st,rare2st).gt.raremin*sdelta.and. & + max(rare1st,rare2st).lt.raremax*sdelta) then + rarecorrector=.true. + if (rare1st.gt.rare2st) then + lambda(2)=s1m + elseif (rare2st.gt.rare1st) then + lambda(2)=s2m + else + lambda(2)=0.5d0*(s1m+s2m) + endif + endif + endif + if (hstarHLL.lt.min(hL,hR)/5.d0) rarecorrector=.false. + endif + +! ## Is this correct 2-wave when rarecorrector == .true. ?? + do mw=1,num_waves + r(1,mw)=1.d0 + r(2,mw)=lambda(mw) + r(3,mw)=(lambda(mw))**2 + enddo + if (.not.rarecorrector) then + lambda(2) = 0.5d0*(lambda(1)+lambda(3)) +! lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1) + r(1,2)=0.d0 + r(2,2)=0.d0 + r(3,2)=1.d0 + endif +! !--------------------------------------------------- + + +! !determine the steady state wave ------------------- + criticaltol = 1.d-6 + deldelh = -delb + deldelphi = -g*0.5d0*(hR+hL)*delb + +! !determine a few quanitites needed for steady state wave if iterated + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + + !iterate to better determine the steady state wave + convergencetol=1.d-6 + do iter=1,maxiter + !determine steady state wave (this will be subtracted from the delta vectors) + if (min(hLstar,hRstar).lt.drytol.and.rarecorrector) then + rarecorrector=.false. + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + lambda(2) = 0.5d0*(lambda(1)+lambda(3)) +! lambda(2) = max(min(0.5d0*(s1m+s2m),sE2),sE1) + r(1,2)=0.d0 + r(2,2)=0.d0 + r(3,2)=1.d0 + endif + + hbar = max(0.5d0*(hLstar+hRstar),0.d0) + s1s2bar = 0.25d0*(uLstar+uRstar)**2 - g*hbar + s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar + +! !find if sonic problem + sonic=.false. + if (abs(s1s2bar).le.criticaltol) sonic=.true. + if (s1s2bar*s1s2tilde.le.criticaltol) sonic=.true. + if (s1s2bar*sE1*sE2.le.criticaltol) sonic = .true. + if (min(abs(sE1),abs(sE2)).lt.criticaltol) sonic=.true. + if (sE1.lt.0.d0.and.s1m.gt.0.d0) sonic = .true. + if (sE2.gt.0.d0.and.s2m.lt.0.d0) sonic = .true. + if ((uL+sqrt(g*hL))*(uR+sqrt(g*hR)).lt.0.d0) sonic=.true. + if ((uL-sqrt(g*hL))*(uR-sqrt(g*hR)).lt.0.d0) sonic=.true. + +! !find jump in h, deldelh + if (sonic) then + deldelh = -delb + else + deldelh = delb*g*hbar/s1s2bar + endif +! !find bounds in case of critical state resonance, or negative states + if (sE1.lt.-criticaltol.and.sE2.gt.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE2) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE1) + elseif (sE1.ge.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE1) + deldelh = max(deldelh,-hL) + elseif (sE2.le.-criticaltol) then + deldelh = min(deldelh,hR) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE2) + endif + +! !find jump in phi, deldelphi + if (sonic) then + deldelphi = -g*hbar*delb + else + deldelphi = -delb*g*hbar*s1s2tilde/s1s2bar + endif +! !find bounds in case of critical state resonance, or negative states + deldelphi=min(deldelphi,g*max(-hLstar*delb,-hRstar*delb)) + deldelphi=max(deldelphi,g*min(-hLstar*delb,-hRstar*delb)) + + del(1)=delh-deldelh + del(2)=delhu + del(3)=delphi-deldelphi + +! !Determine determinant of eigenvector matrix======== + det1=r(1,1)*(r(2,2)*r(3,3)-r(2,3)*r(3,2)) + det2=r(1,2)*(r(2,1)*r(3,3)-r(2,3)*r(3,1)) + det3=r(1,3)*(r(2,1)*r(3,2)-r(2,2)*r(3,1)) + determinant=det1-det2+det3 + +! !solve for beta(k) using Cramers Rule================= + do k=1,3 + do mw=1,3 + A(1,mw)=r(1,mw) + A(2,mw)=r(2,mw) + A(3,mw)=r(3,mw) + enddo + A(1,k)=del(1) + A(2,k)=del(2) + A(3,k)=del(3) + det1=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) + det2=A(1,2)*(A(2,1)*A(3,3)-A(2,3)*A(3,1)) + det3=A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) + beta(k)=(det1-det2+det3)/determinant + enddo + + !exit if things aren't changing + if (abs(del(1)**2+del(3)**2-delnorm).lt.convergencetol) exit + delnorm = del(1)**2+del(3)**2 + !find new states qLstar and qRstar on either side of interface + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + huLstar=uLstar*hLstar + huRstar=uRstar*hRstar + do mw=1,num_waves + if (lambda(mw).lt.0.d0) then + hLstar= hLstar + beta(mw)*r(1,mw) + huLstar= huLstar + beta(mw)*r(2,mw) + endif + enddo + do mw=num_waves,1,-1 + if (lambda(mw).gt.0.d0) then + hRstar= hRstar - beta(mw)*r(1,mw) + huRstar= huRstar - beta(mw)*r(2,mw) + endif + enddo + + if (hLstar.gt.drytol) then + uLstar=huLstar/hLstar + else + hLstar=max(hLstar,0.d0) + uLstar=0.d0 + endif + if (hRstar.gt.drytol) then + uRstar=huRstar/hRstar + else + hRstar=max(hRstar,0.d0) + uRstar=0.d0 + endif + + enddo ! end iteration on Riemann problem + + do mw=1,num_waves + sw(mw)=lambda(mw) + fw(1,mw)=beta(mw)*r(2,mw) + fw(2,mw)=beta(mw)*r(3,mw) + fw(3,mw)=beta(mw)*r(2,mw) + enddo + !find transverse components (ie huv jumps). + fw(3,1)=fw(3,1)*vL + fw(3,3)=fw(3,3)*vR + fw(3,2)= hR*uR*vR - hL*uL*vL - fw(3,1)- fw(3,3) + + return + + end !subroutine riemann_aug_JCP------------------------------------------------- + + +!----------------------------------------------------------------------- + subroutine riemann_ssqfwave(maxiter,num_eqn,num_waves,hL,hR,huL,huR, & + hvL,hvR,bL,bR,uL,uR,vL,vR,phiL,phiR,sE1,sE2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! steady state wave is subtracted from delta [q,f]^T before decomposition + + implicit none + + !input + integer num_eqn,num_waves,maxiter + + real(kind=8) hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,sE1,sE2 + real(kind=8) vL,vR,hvL,hvR + real(kind=8) drytol,g + + !local + integer iter + + logical sonic + + real(kind=8) delh,delhu,delphi,delb,delhdecomp,delphidecomp + real(kind=8) s1s2bar,s1s2tilde,hbar,hLstar,hRstar,hustar + real(kind=8) uRstar,uLstar,hstarHLL + real(kind=8) deldelh,deldelphi + real(kind=8) alpha1,alpha2,beta1,beta2,delalpha1,delalpha2 + real(kind=8) criticaltol,convergencetol + real(kind=8) sL,sR + real(kind=8) uhat,chat,sRoe1,sRoe2 + + real(kind=8) sw(num_waves) + real(kind=8) fw(num_eqn,num_waves) + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + + convergencetol= 1.d-16 + criticaltol = 1.d-99 + + deldelh = -delb + deldelphi = -g*0.5d0*(hR+hL)*delb + +! !if no source term, skip determining steady state wave + if (abs(delb).gt.0.d0) then +! + !determine a few quanitites needed for steady state wave if iterated + hLstar=hL + hRstar=hR + uLstar=uL + uRstar=uR + hstarHLL = max((huL-huR+sE2*hR-sE1*hL)/(sE2-sE1),0.d0) ! middle state in an HLL solve + + alpha1=0.d0 + alpha2=0.d0 + +! !iterate to better determine Riemann problem + do iter=1,maxiter + + !determine steady state wave (this will be subtracted from the delta vectors) + hbar = max(0.5d0*(hLstar+hRstar),0.d0) + s1s2bar = 0.25d0*(uLstar+uRstar)**2 - g*hbar + s1s2tilde= max(0.d0,uLstar*uRstar) - g*hbar + + +! !find if sonic problem + sonic=.false. + if (abs(s1s2bar).le.criticaltol) sonic=.true. + if (s1s2bar*s1s2tilde.le.criticaltol) sonic=.true. + if (s1s2bar*sE1*sE2.le.criticaltol) sonic = .true. + if (min(abs(sE1),abs(sE2)).lt.criticaltol) sonic=.true. + +! !find jump in h, deldelh + if (sonic) then + deldelh = -delb + else + deldelh = delb*g*hbar/s1s2bar + endif +! !bounds in case of critical state resonance, or negative states + if (sE1.lt.-criticaltol.and.sE2.gt.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE2) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE1) + elseif (sE1.ge.criticaltol) then + deldelh = min(deldelh,hstarHLL*(sE2-sE1)/sE1) + deldelh = max(deldelh,-hL) + elseif (sE2.le.-criticaltol) then + deldelh = min(deldelh,hR) + deldelh = max(deldelh,hstarHLL*(sE2-sE1)/sE2) + endif + +! !find jump in phi, deldelphi + if (sonic) then + deldelphi = -g*hbar*delb + else + deldelphi = -delb*g*hbar*s1s2tilde/s1s2bar + endif +! !bounds in case of critical state resonance, or negative states + deldelphi=min(deldelphi,g*max(-hLstar*delb,-hRstar*delb)) + deldelphi=max(deldelphi,g*min(-hLstar*delb,-hRstar*delb)) + +!---------determine fwaves ------------------------------------------ + +! !first decomposition + delhdecomp = delh-deldelh + delalpha1 = (sE2*delhdecomp - delhu)/(sE2-sE1)-alpha1 + alpha1 = alpha1 + delalpha1 + delalpha2 = (delhu - sE1*delhdecomp)/(sE2-sE1)-alpha2 + alpha2 = alpha2 + delalpha2 + + !second decomposition + delphidecomp = delphi - deldelphi + beta1 = (sE2*delhu - delphidecomp)/(sE2-sE1) + beta2 = (delphidecomp - sE1*delhu)/(sE2-sE1) + + if ((delalpha2**2+delalpha1**2).lt.convergencetol**2) then + exit + endif +! + if (sE2.gt.0.d0.and.sE1.lt.0.d0) then + hLstar=hL+alpha1 + hRstar=hR-alpha2 +! hustar=huL+alpha1*sE1 + hustar = huL + beta1 + elseif (sE1.ge.0.d0) then + hLstar=hL + hustar=huL + hRstar=hR - alpha1 - alpha2 + elseif (sE2.le.0.d0) then + hRstar=hR + hustar=huR + hLstar=hL + alpha1 + alpha2 + endif +! + if (hLstar.gt.drytol) then + uLstar=hustar/hLstar + else + hLstar=max(hLstar,0.d0) + uLstar=0.d0 + endif +! + if (hRstar.gt.drytol) then + uRstar=hustar/hRstar + else + hRstar=max(hRstar,0.d0) + uRstar=0.d0 + endif + + enddo + endif + + delhdecomp = delh - deldelh + delphidecomp = delphi - deldelphi + + !first decomposition + alpha1 = (sE2*delhdecomp - delhu)/(sE2-sE1) + alpha2 = (delhu - sE1*delhdecomp)/(sE2-sE1) + + !second decomposition + beta1 = (sE2*delhu - delphidecomp)/(sE2-sE1) + beta2 = (delphidecomp - sE1*delhu)/(sE2-sE1) + + ! 1st nonlinear wave + fw(1,1) = alpha1*sE1 + fw(2,1) = beta1*sE1 + fw(3,1) = fw(1,1)*vL + ! 2nd nonlinear wave + fw(1,3) = alpha2*sE2 + fw(2,3) = beta2*sE2 + fw(3,3) = fw(1,3)*vR + ! advection of transverse wave + fw(1,2) = 0.d0 + fw(2,2) = 0.d0 + fw(3,2) = hR*uR*vR - hL*uL*vL -fw(3,1)-fw(3,3) + !speeds + sw(1)=sE1 + sw(2)=0.5d0*(sE1+sE2) + sw(3)=sE2 + + return + + end subroutine !------------------------------------------------- + + +!----------------------------------------------------------------------- + subroutine riemann_fwave(num_eqn,num_waves,hL,hR,huL,huR,hvL,hvR, & + bL,bR,uL,uR,vL,vR,phiL,phiR,s1,s2,drytol,g,sw,fw) + + ! solve shallow water equations given single left and right states + ! solution has two waves. + ! flux - source is decomposed. + + implicit none + + !input + integer num_eqn,num_waves + + real(kind=8) hL,hR,huL,huR,bL,bR,uL,uR,phiL,phiR,s1,s2 + real(kind=8) hvL,hvR,vL,vR + real(kind=8) drytol,g + + real(kind=8) sw(num_waves) + real(kind=8) fw(num_eqn,num_waves) + + !local + real(kind=8) delh,delhu,delphi,delb,delhdecomp,delphidecomp + real(kind=8) deldelh,deldelphi + real(kind=8) beta1,beta2 + + + !determine del vectors + delh = hR-hL + delhu = huR-huL + delphi = phiR-phiL + delb = bR-bL + + deldelphi = -g*0.5d0*(hR+hL)*delb + delphidecomp = delphi - deldelphi + + !flux decomposition + beta1 = (s2*delhu - delphidecomp)/(s2-s1) + beta2 = (delphidecomp - s1*delhu)/(s2-s1) + + sw(1)=s1 + sw(2)=0.5d0*(s1+s2) + sw(3)=s2 + ! 1st nonlinear wave + fw(1,1) = beta1 + fw(2,1) = beta1*s1 + fw(3,1) = beta1*vL + ! 2nd nonlinear wave + fw(1,3) = beta2 + fw(2,3) = beta2*s2 + fw(3,3) = beta2*vR + ! advection of transverse wave + fw(1,2) = 0.d0 + fw(2,2) = 0.d0 + fw(3,2) = hR*uR*vR - hL*uL*vL -fw(3,1)-fw(3,3) + return + + end !subroutine ------------------------------------------------- + + + + + +!============================================================================= + subroutine riemanntype(hL,hR,uL,uR,hm,s1m,s2m,rare1,rare2, & + maxiter,drytol,g) + + !determine the Riemann structure (wave-type in each family) + + + implicit none + + !input + real(kind=8) hL,hR,uL,uR,drytol,g + integer maxiter + + !output + real(kind=8) s1m,s2m + logical rare1,rare2 + + !local + real(kind=8) hm,u1m,u2m,um,delu + real(kind=8) h_max,h_min,h0,F_max,F_min,dfdh,F0,slope,gL,gR + integer iter + + + +! !Test for Riemann structure + + h_min=min(hR,hL) + h_max=max(hR,hL) + delu=uR-uL + + if (h_min.le.drytol) then + hm=0.d0 + um=0.d0 + s1m=uR+uL-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hL) + s2m=uR+uL-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hL) + if (hL.le.0.d0) then + rare2=.true. + rare1=.false. + else + rare1=.true. + rare2=.false. + endif + + else + F_min= delu+2.d0*(sqrt(g*h_min)-sqrt(g*h_max)) + F_max= delu + (h_max-h_min)*(sqrt(.5d0*g*(h_max+h_min)/(h_max*h_min))) + + if (F_min.gt.0.d0) then !2-rarefactions + + hm=(1.d0/(16.d0*g))* & + max(0.d0,-delu+2.d0*(sqrt(g*hL)+sqrt(g*hR)))**2 + um=sign(1.d0,hm)*(uL+2.d0*(sqrt(g*hL)-sqrt(g*hm))) + + s1m=uL+2.d0*sqrt(g*hL)-3.d0*sqrt(g*hm) + s2m=uR-2.d0*sqrt(g*hR)+3.d0*sqrt(g*hm) + + rare1=.true. + rare2=.true. + + elseif (F_max.le.0.d0) then !2 shocks + +! !root finding using a Newton iteration on sqrt(h)=== + h0=h_max + do iter=1,maxiter + gL=sqrt(.5d0*g*(1/h0 + 1/hL)) + gR=sqrt(.5d0*g*(1/h0 + 1/hR)) + F0=delu+(h0-hL)*gL + (h0-hR)*gR + dfdh=gL-g*(h0-hL)/(4.d0*(h0**2)*gL)+gR-g*(h0-hR)/(4.d0*(h0**2)*gR) + slope=2.d0*sqrt(h0)*dfdh + h0=(sqrt(h0)-F0/slope)**2 + enddo + hm=h0 + u1m=uL-(hm-hL)*sqrt((.5d0*g)*(1/hm + 1/hL)) + u2m=uR+(hm-hR)*sqrt((.5d0*g)*(1/hm + 1/hR)) + um=.5d0*(u1m+u2m) + + s1m=u1m-sqrt(g*hm) + s2m=u2m+sqrt(g*hm) + rare1=.false. + rare2=.false. + + else !one shock one rarefaction + h0=h_min + + do iter=1,maxiter + F0=delu + 2.d0*(sqrt(g*h0)-sqrt(g*h_max)) & + + (h0-h_min)*sqrt(.5d0*g*(1/h0+1/h_min)) + slope=(F_max-F0)/(h_max-h_min) + h0=h0-F0/slope + enddo + + hm=h0 + if (hL.gt.hR) then + um=uL+2.d0*sqrt(g*hL)-2.d0*sqrt(g*hm) + s1m=uL+2.d0*sqrt(g*hL)-3.d0*sqrt(g*hm) + s2m=uL+2.d0*sqrt(g*hL)-sqrt(g*hm) + rare1=.true. + rare2=.false. + else + s2m=uR-2.d0*sqrt(g*hR)+3.d0*sqrt(g*hm) + s1m=uR-2.d0*sqrt(g*hR)+sqrt(g*hm) + um=uR-2.d0*sqrt(g*hR)+2.d0*sqrt(g*hm) + rare2=.true. + rare1=.false. + endif + endif + endif + + return + + end ! subroutine riemanntype---------------------------------------------------------------- diff --git a/src/rpt2_shallow_bathymetry_fwave.f90 b/src/rpt2_shallow_bathymetry_fwave.f90 new file mode 100755 index 00000000..794dd6e7 --- /dev/null +++ b/src/rpt2_shallow_bathymetry_fwave.f90 @@ -0,0 +1,119 @@ +subroutine rpt2(ixy, imp, maxm, num_eqn, num_waves, num_aux, num_ghost, & + num_cells, ql, qr, aux1, aux2, aux3, asdq, bmasdq, bpasdq) +! Transverse solver for shallow water with bathymetry using an fwave solver +! Note that this solver does NOT handle dry-states +! +! (h)_t + (h*u)_x = 0 +! (hu)_t + (h*u^2 + 1/2*grav*h^2)_x + (h*u*v)_y = -grav*h*(b)_x +! (hv)_t + (h*u*v)_x + (h*v^2 + 1/2*grav*h^2)_y = -grav*h*(b)_y +! +! using the f-wave algorithm and Roe's approximate Riemann solver. +! +! waves: 3 +! equations: 3 +! +! Conserved quantities: +! 1 depth +! 2 x_momentum +! 3 y_momentum +! +! Auxiliary fields: +! 1 bathymetry +! +! The gravitational constant grav should be in the common block cparam. +! +! See http://www.clawpack.org/riemann.html for a detailed explanation +! of the Riemann solver API. + + implicit none + + ! Input + integer, intent(in) :: ixy, imp, maxm, num_eqn, num_waves, num_aux + integer, intent(in) :: num_ghost, num_cells + real(kind=8), intent(in) :: ql(num_eqn, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(in) :: qr(num_eqn, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(in) :: aux1(num_aux, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(in) :: aux2(num_aux, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(in) :: aux3(num_aux, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(in) :: asdq(num_eqn, 1 - num_ghost:maxm + num_ghost) + + ! Output + real(kind=8), intent(out) :: bmasdq(num_eqn, 1 - num_ghost:maxm + num_ghost) + real(kind=8), intent(out) :: bpasdq(num_eqn, 1 - num_ghost:maxm + num_ghost) + + ! Local + integer :: i, k, normal_index, transverse_index + real(kind=8) :: h_l, h_r, u_l, u_r, v_l, v_r, h_bar, h_root, u_hat, v_hat + real(kind=8) :: dry_state_l, dry_state_r, s(3), beta(3), R(3, 3) + + ! Parameters + real(kind=8) :: grav, dry_tolerance, sea_level + common /cparam/ grav, dry_tolerance, sea_level + + ! Deterimine direction of solve + if (ixy == 1) then + normal_index = 2 + transverse_index = 3 + else + normal_index = 3 + transverse_index = 2 + end if + + ! Primary loop over cell edges + do i = 2 - num_ghost, num_cells + num_ghost + + ! Check for dry states - need merge here to convert to float + dry_state_l = merge(0.d0, 1.d0, qr(1, i - 1) < dry_tolerance) + dry_state_r = merge(0.d0, 1.d0, ql(1, i) < dry_tolerance) + + ! Extract state variables + h_l = qr(1, i - 1) * dry_state_l + u_l = qr(normal_index, i - 1) / qr(1, i - 1) * dry_state_l + v_l = qr(transverse_index, i - 1) / qr(1, i - 1) * dry_state_l + + h_r = ql(1, i) * dry_state_r + u_r = ql(normal_index, i) / ql(1, i) * dry_state_r + v_r = ql(transverse_index, i) / ql(1, i) * dry_state_r + + ! Determine speeds for the Jacobian + h_bar = 0.5d0 * (h_r + h_l) + h_root = sqrt(h_r) + sqrt(h_l) + v_hat = (v_r * sqrt(h_r)) / h_root + (v_l * sqrt(h_l)) / h_root + u_hat = (u_r * sqrt(h_r)) / h_root + (u_l * sqrt(h_l)) / h_root + + s(1) = min(v_hat - sqrt(grav * h_bar), v_l - sqrt(grav * h_l)) + s(3) = max(v_hat + sqrt(grav * h_bar), v_r + sqrt(grav * h_r)) + s(2) = 0.5d0 * (s(1) + s(2)) + + ! Determine asdq decompisition + beta(1) = s(3) * asdq(1, i) / (s(3) - s(1)) & + - asdq(transverse_index, i) / (s(3) - s(1)) + beta(2) = -s(2) * asdq(1, i) + asdq(normal_index, i) + beta(3) = asdq(transverse_index, i) / (s(3) - s(1)) & + - (s(1) * asdq(1, i) / (s(3) - s(1))) + + ! Eigenvectors + R(1, :) = [1.d0, 0.d0, 1.d0] + R(2, :) = [s(2), 1.d0, s(2)] + R(3, :) = [s(1), 0.d0, s(3)] + + ! Compute fluctuations + do k = 1, num_waves + if (s(k) < 0.d0) then + bmasdq(1, i) = bmasdq(1, i) + s(k) * beta(k) * R(1, k) + bmasdq(normal_index, i) = bmasdq(normal_index, i) & + + s(k) * beta(k) * R(2, k) + bmasdq(transverse_index, i) = bmasdq(transverse_index, i) & + + s(k) * beta(k) * R(3, k) + else + bpasdq(1, i) = bpasdq(1, i) + s(k) * beta(k) * R(1, k) + bpasdq(normal_index, i) = bpasdq(normal_index, i) & + + s(k) * beta(k) * R(2, k) + bpasdq(transverse_index, i) = bpasdq(transverse_index, i) & + + s(k) * beta(k) * R(3, k) + end if + end do + + end do ! End of primary loop + +end subroutine rpt2 \ No newline at end of file diff --git a/src/rpt2_sw_aug.f90 b/src/rpt2_sw_aug.f90 new file mode 100755 index 00000000..343be65e --- /dev/null +++ b/src/rpt2_sw_aug.f90 @@ -0,0 +1,178 @@ +! ===================================================== +subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq) +! ===================================================== + + implicit none +! +! # Riemann solver in the transverse direction using an einfeldt +! Jacobian. + + double precision :: grav, g + double precision, parameter :: tol = 1.e-14 + common /cparam/ grav + + integer ixy,maxm,meqn,maux,mwaves,mbc,mx,imp + + double precision ql(meqn,1-mbc:maxm+mbc) + double precision qr(meqn,1-mbc:maxm+mbc) + double precision asdq(meqn,1-mbc:maxm+mbc) + double precision bmasdq(meqn,1-mbc:maxm+mbc) + double precision bpasdq(meqn,1-mbc:maxm+mbc) + double precision aux1(maux,1-mbc:maxm+mbc) + double precision aux2(maux,1-mbc:maxm+mbc) + double precision aux3(maux,1-mbc:maxm+mbc) + + double precision s(3) + double precision r(3,3) + double precision beta(3) + double precision abs_tol + double precision hl,hr,hul,hur,hvl,hvr,vl,vr,ul,ur,bl,br + double precision uhat,vhat,hhat,roe1,roe3,s1,s2,s3,s1l,s3r + double precision delf1,delf2,delf3,dxdcd,dxdcu + double precision dxdcm,dxdcp,topo1,topo3,eta + + integer i,m,mw,mu,mv + + abs_tol=tol + g = grav + + if (ixy.eq.1) then + mu = 2 + mv = 3 + else + mu = 3 + mv = 2 + endif + + do i=2-mbc,mx+mbc + + hl=qr(1,i-1) + hr=ql(1,i) + hul=qr(mu,i-1) + hur=ql(mu,i) + hvl=qr(mv,i-1) + hvr=ql(mv,i) + +!===========determine velocity from momentum=========================== + if (hl.lt.abs_tol) then + hl=0.d0 + ul=0.d0 + vl=0.d0 + else + ul=hul/hl + vl=hvl/hl + endif + + if (hr.lt.abs_tol) then + hr=0.d0 + ur=0.d0 + vr=0.d0 + else + ur=hur/hr + vr=hvr/hr + endif + + do mw=1,mwaves + s(mw)=0.d0 + beta(mw)=0.d0 + do m=1,meqn + r(m,mw)=0.d0 + enddo + enddo + dxdcp = 1.d0 + dxdcm = 1.d0 + + if (hl <= tol .and. hr <= tol) go to 90 + + !check and see if cell that transverse waves are going in is high and dry + if (imp.eq.1) then + eta = qr(1,i-1) + aux2(1,i-1) + topo1 = aux1(1,i-1) + topo3 = aux3(1,i-1) +! s1 = vl-sqrt(g*hl) +! s3 = vl+sqrt(g*hl) +! s2 = 0.5d0*(s1+s3) + else + eta = ql(1,i) + aux2(1,i) + topo1 = aux1(1,i) + topo3 = aux3(1,i) +! s1 = vr-sqrt(g*hr) +! s3 = vr+sqrt(g*hr) +! s2 = 0.5d0*(s1+s3) + endif + if (eta.lt.max(topo1,topo3)) go to 90 + + +!=====Determine some speeds necessary for the Jacobian================= + vhat=(vr*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + & + (vl*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) + + uhat=(ur*dsqrt(hr))/(dsqrt(hr)+dsqrt(hl)) + & + (ul*dsqrt(hl))/(dsqrt(hr)+dsqrt(hl)) + hhat=(hr+hl)/2.d0 + + roe1=vhat-dsqrt(g*hhat) + roe3=vhat+dsqrt(g*hhat) + + s1l=vl-dsqrt(g*hl) + s3r=vr+dsqrt(g*hr) + + s1=dmin1(roe1,s1l) + s3=dmax1(roe3,s3r) + + s2=0.5d0*(s1+s3) + + s(1)=s1 + s(2)=s2 + s(3)=s3 +!=======================Determine asdq decomposition (beta)============ + delf1=asdq(1,i) + delf2=asdq(mu,i) + delf3=asdq(mv, i) + + beta(1) = (s3*delf1/(s3-s1))-(delf3/(s3-s1)) + beta(2) = -s2*delf1 + delf2 + beta(3) = (delf3/(s3-s1))-(s1*delf1/(s3-s1)) +!======================End ================================================= + +!=====================Set-up eigenvectors=================================== + r(1,1) = 1.d0 + r(2,1) = s2 + r(3,1) = s1 + + r(1,2) = 0.d0 + r(2,2) = 1.d0 + r(3,2) = 0.d0 + + r(1,3) = 1.d0 + r(2,3) = s2 + r(3,3) = s3 +!============================================================================ +90 continue +!============= compute fluctuations========================================== + + bmasdq(1,i)=0.0d0 + bpasdq(1,i)=0.0d0 + bmasdq(2,i)=0.0d0 + bpasdq(2,i)=0.0d0 + bmasdq(3,i)=0.0d0 + bpasdq(3,i)=0.0d0 + do mw=1,3 + if (s(mw).lt.0.d0) then + bmasdq(1,i) =bmasdq(1,i) + dxdcm*s(mw)*beta(mw)*r(1,mw) + bmasdq(mu,i)=bmasdq(mu,i)+ dxdcm*s(mw)*beta(mw)*r(2,mw) + bmasdq(mv,i)=bmasdq(mv,i)+ dxdcm*s(mw)*beta(mw)*r(3,mw) + elseif (s(mw).gt.0.d0) then + bpasdq(1,i) =bpasdq(1,i) + dxdcp*s(mw)*beta(mw)*r(1,mw) + bpasdq(mu,i)=bpasdq(mu,i)+ dxdcp*s(mw)*beta(mw)*r(2,mw) + bpasdq(mv,i)=bpasdq(mv,i)+ dxdcp*s(mw)*beta(mw)*r(3,mw) + endif + enddo +!======================================================================== + enddo +! + +! + + return + end diff --git a/src/setup.py b/src/setup.py index 075fbaef..04812fb2 100644 --- a/src/setup.py +++ b/src/setup.py @@ -19,6 +19,7 @@ 'nonlinear_elasticity_fwave', 'reactive_euler_with_efix', 'shallow_roe_with_efix', + 'shallow_bathymetry_fwave', 'shallow_roe_tracer'] two_d_ptwise_riemann = ['acoustics'] @@ -30,6 +31,8 @@ 'euler_5wave', 'psystem', 'shallow_roe_with_efix', + 'shallow_bathymetry_fwave', + 'sw_aug', 'shallow_sphere', 'vc_acoustics', 'vc_advection', diff --git a/src/shallow_1D_py.py b/src/shallow_1D_py.py index bc010cdd..1cc05212 100644 --- a/src/shallow_1D_py.py +++ b/src/shallow_1D_py.py @@ -2,39 +2,29 @@ # encoding: utf-8 r""" Riemann solvers for the shallow water equations. - + The available solvers are: * Roe - Use Roe averages to caluclate the solution to the Riemann problem * HLL - Use a HLL solver - * Exact - Use a newton iteration to calculate the exact solution to the + * Exact - Use a newton iteration to calculate the exact solution to the Riemann problem -.. math:: +.. math:: q_t + f(q)_x = 0 - -where -.. math:: +where + +.. math:: q(x,t) = \left [ \begin{array}{c} h \\ h u \end{array} \right ], - -the flux function is -.. math:: +the flux function is + +.. math:: f(q) = \left [ \begin{array}{c} h u \\ hu^2 + 1/2 g h^2 \end{array}\right ]. and :math:`h` is the water column height, :math:`u` the velocity and :math:`g` is the gravitational acceleration. - -:Authors: - Kyle T. Mandli (2009-02-05): Initial version """ -# ============================================================================ -# Copyright (C) 2009 Kyle T. Mandli -# -# Distributed under the terms of the Berkeley Software Distribution (BSD) -# license -# http://www.opensource.org/licenses/ -# ============================================================================ from __future__ import absolute_import import numpy as np @@ -43,33 +33,34 @@ num_eqn = 2 num_waves = 2 -def shallow_roe_1D(q_l,q_r,aux_l,aux_r,problem_data): + +def shallow_roe_1D(q_l, q_r, aux_l, aux_r, problem_data): r""" Roe shallow water solver in 1d:: - + ubar = (sqrt(u_l) + sqrt(u_r)) / (sqrt(h_l) + sqrt(h_r)) cbar = sqrt( 0.5 * g * (h_l + h_r)) - + W_1 = | 1 | s_1 = ubar - cbar | ubar - cbar | - + W_2 = | 1 | s_1 = ubar + cbar | ubar + cbar | - + a1 = 0.5 * ( - delta_hu + (ubar + cbar) * delta_h ) / cbar a2 = 0.5 * ( delta_hu - (ubar - cbar) * delta_h ) / cbar - + *problem_data* should contain: - *g* - (float) Gravitational constant - - *efix* - (bool) Boolean as to whether a entropy fix should be used, if + - *efix* - (bool) Boolean as to whether a entropy fix should be used, if not present, false is assumed - + :Version: 1.0 (2009-02-05) """ - + # Array shapes num_rp = q_l.shape[1] - + # Output arrays wave = np.empty( (num_eqn, num_waves, num_rp) ) s = np.zeros( (num_waves, num_rp) ) @@ -78,23 +69,23 @@ def shallow_roe_1D(q_l,q_r,aux_l,aux_r,problem_data): # Compute roe-averaged quantities ubar = ( (q_l[1,:]/np.sqrt(q_l[0,:]) + q_r[1,:]/np.sqrt(q_r[0,:])) / - (np.sqrt(q_l[0,:]) + np.sqrt(q_r[0,:])) ) + (np.sqrt(q_l[0,:]) + np.sqrt(q_r[0,:])) ) cbar = np.sqrt(0.5 * problem_data['grav'] * (q_l[0,:] + q_r[0,:])) - - # Compute Flux structure + + # Compute Flux structure delta = q_r - q_l a1 = 0.5 * (-delta[1,:] + (ubar + cbar) * delta[0,:]) / cbar a2 = 0.5 * ( delta[1,:] - (ubar - cbar) * delta[0,:]) / cbar - + # Compute each family of waves wave[0,0,:] = a1 wave[1,0,:] = a1 * (ubar - cbar) s[0,:] = ubar - cbar - + wave[0,1,:] = a2 wave[1,1,:] = a2 * (ubar + cbar) s[1,:] = ubar + cbar - + if problem_data['efix']: raise NotImplementedError("Entropy fix has not been implemented.") else: @@ -104,29 +95,29 @@ def shallow_roe_1D(q_l,q_r,aux_l,aux_r,problem_data): s_index[0,:] = s[mw,:] amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:] apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:] - + return wave, s, amdq, apdq - + def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data): r""" HLL shallow water solver :: - - + + W_1 = Q_hat - Q_l s_1 = min(u_l-c_l,u_l+c_l,lambda_roe_1,lambda_roe_2) W_2 = Q_r - Q_hat s_2 = max(u_r-c_r,u_r+c_r,lambda_roe_1,lambda_roe_2) - + Q_hat = ( f(q_r) - f(q_l) - s_2 * q_r + s_1 * q_l ) / (s_1 - s_2) - + *problem_data* should contain: - *g* - (float) Gravitational constant - + :Version: 1.0 (2009-02-05) """ # Array shapes num_rp = q_l.shape[1] num_eqn = 2 num_waves = 2 - + # Output arrays wave = np.empty( (num_eqn, num_waves, num_rp) ) s = np.empty( (num_waves, num_rp) ) @@ -155,7 +146,7 @@ def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data): # Compute middle state q_hat = np.empty((2,num_rp)) - q_hat[0,:] = ((q_r[1,:] - q_l[1,:] - s[1,:] * q_r[0,:] + q_hat[0,:] = ((q_r[1,:] - q_l[1,:] - s[1,:] * q_r[0,:] + s[0,:] * q_l[0,:]) / (s[0,:] - s[1,:])) q_hat[1,:] = ((q_r[1,:]**2/q_r[0,:] + 0.5 * problem_data['grav'] * q_r[0,:]**2 - (q_l[1,:]**2/q_l[0,:] + 0.5 * problem_data['grav'] * q_l[0,:]**2) @@ -164,7 +155,7 @@ def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data): # Compute each family of waves wave[:,0,:] = q_hat - q_l wave[:,1,:] = q_r - q_hat - + # Compute variations s_index = np.zeros((2,num_rp)) for m in range(num_eqn): @@ -172,7 +163,7 @@ def shallow_hll_1D(q_l,q_r,aux_l,aux_r,problem_data): s_index[0,:] = s[mw,:] amdq[m,:] += np.min(s_index,axis=0) * wave[m,mw,:] apdq[m,:] += np.max(s_index,axis=0) * wave[m,mw,:] - + return wave, s, amdq, apdq @@ -180,21 +171,26 @@ def shallow_fwave_1d(q_l, q_r, aux_l, aux_r, problem_data): r"""Shallow water Riemann solver using fwaves Also includes support for bathymetry but be wary if you think you might have - dry states as this has not been tested. - + dry states as this has not been tested. + *problem_data* should contain: - *grav* - (float) Gravitational constant + - *dry_tolerance* - (float) Set velocities to zero if h is below this + tolerance. - *sea_level* - (float) Datum from which the dry-state is calculated. - + :Version: 1.0 (2014-09-05) + :Version: 2.0 (2017-03-07) """ g = problem_data['grav'] + dry_tolerance = problem_data['dry_tolerance'] + sea_level = problem_data['sea_level'] num_rp = q_l.shape[1] num_eqn = 2 num_waves = 2 - + # Output arrays fwave = np.empty( (num_eqn, num_waves, num_rp) ) s = np.empty( (num_waves, num_rp) ) @@ -202,42 +198,51 @@ def shallow_fwave_1d(q_l, q_r, aux_l, aux_r, problem_data): apdq = np.zeros( (num_eqn, num_rp) ) # Extract state - u_l = np.where(q_l[0,:] - problem_data['sea_level'] > 1e-3, - q_l[1,:] / q_l[0,:], 0.0) - u_r = np.where(q_r[0,:] - problem_data['sea_level'] > 1e-3, - q_r[1,:] / q_r[0,:], 0.0) - phi_l = q_l[0,:] * u_l**2 + 0.5 * g * q_l[0,:]**2 - phi_r = q_r[0,:] * u_r**2 + 0.5 * g * q_r[0,:]**2 + u_l = np.where(q_l[0, :] > dry_tolerance, + q_l[1, :] / q_l[0, :], 0.0) + u_r = np.where(q_r[0, :] > dry_tolerance, + q_r[1, :] / q_r[0, :], 0.0) + phi_l = q_l[0, :] * u_l**2 + 0.5 * g * q_l[0, :]**2 + phi_r = q_r[0, :] * u_r**2 + 0.5 * g * q_r[0, :]**2 + h_bar = 0.5 * (q_l[0, :] + q_r[0, :]) # Speeds - s[0,:] = u_l - np.sqrt(g * q_l[0,:]) - s[1,:] = u_r + np.sqrt(g * q_r[0,:]) + u_hat = (np.sqrt(g * q_l[0, :]) * u_l + np.sqrt(g * q_r[0, :]) * u_r) \ + / (np.sqrt(g * q_l[0, :]) + np.sqrt(g * q_r[0, :])) + c_hat = np.sqrt(g * h_bar) + s[0, :] = np.amin(np.vstack((u_l - np.sqrt(g * q_l[0, :]), + u_hat - c_hat)), axis=0) + s[1, :] = np.amax(np.vstack((u_r + np.sqrt(g * q_r[0, :]), + u_hat + c_hat)), axis=0) - delta1 = q_r[1,:] - q_l[1,:] - delta2 = phi_r - phi_l + g * 0.5 * (q_r[0,:] + q_l[0,:]) * (aux_r[0,:] - aux_l[0,:]) + delta1 = q_r[1, :] - q_l[1, :] + delta2 = phi_r - phi_l + g * h_bar * (aux_r[0, :] - aux_l[0, :]) - beta1 = (s[1,:] * delta1 - delta2) / (s[1,:] - s[0,:]) - beta2 = (delta2 - s[0,:] * delta1) / (s[1,:] - s[0,:]) + beta1 = (s[1, :] * delta1 - delta2) / (s[1, :] - s[0, :]) + beta2 = (delta2 - s[0, :] * delta1) / (s[1, :] - s[0, :]) - fwave[0,0,:] = beta1 - fwave[1,0,:] = beta1 * s[0,:] - fwave[0,1,:] = beta2 - fwave[1,1,:] = beta2 * s[1,:] + fwave[0, 0, :] = beta1 + fwave[1, 0, :] = beta1 * s[0, :] + fwave[0, 1, :] = beta2 + fwave[1, 1, :] = beta2 * s[1, :] for m in range(num_eqn): for mw in range(num_waves): - amdq[m,:] += (s[mw,:] < 0.0) * fwave[m,mw,:] - apdq[m,:] += (s[mw,:] >= 0.0) * fwave[m,mw,:] + amdq[m, :] += (s[mw, :] < 0.0) * fwave[m, mw, :] + apdq[m, :] += (s[mw, :] > 0.0) * fwave[m, mw, :] + + amdq[m, :] += (s[mw, :] == 0.0) * fwave[m, mw, :] * 0.5 + apdq[m, :] += (s[mw, :] == 0.0) * fwave[m, mw, :] * 0.5 return fwave, s, amdq, apdq - -def shallow_exact_1D(q_l,q_r,aux_l,aux_r,problem_data): + +def shallow_exact_1D(q_l, q_r, aux_l, aux_r, problem_data): r""" Exact shallow water Riemann solver - + .. warning:: This solver has not been implemented. - + """ raise NotImplementedError("The exact swe solver has not been implemented.") diff --git a/src/static.py b/src/static.py index b31a2642..4b96e721 100644 --- a/src/static.py +++ b/src/static.py @@ -20,6 +20,8 @@ 'shallow_roe_tracer_1D' : 3, 'shallow_roe_1D' : 2, 'shallow_hll_1D' : 2, + 'shallow_bathymetry_fwave_1D' : 2, + 'shallow_1D_py' : 2, 'vc_advection_1D' : 1, 'layered_shallow_water_1D' : 4, 'acoustics_2D' : 3, @@ -31,6 +33,8 @@ 'kpp_2D' : 1, 'psystem_2D' : 3, 'shallow_roe_with_efix_2D' : 3, + 'shallow_bathymetry_fwave_2D' : 3, + 'sw_aug_2D' : 3, 'shallow_sphere_2D' : 4, 'vc_advection_2D' : 1, 'vc_acoustics_2D' : 3, @@ -61,6 +65,8 @@ 'shallow_roe_tracer_1D' : 3, 'shallow_roe_1D' : 2, 'shallow_hll_1D' : 2, + 'shallow_bathymetry_fwave_1D' : 2, + 'shallow_1D_py' : 2, 'vc_advection_1D' : 1, 'layered_shallow_water_1D' : 4, 'acoustics_2D' : 2, @@ -72,6 +78,8 @@ 'kpp_2D' : 1, 'psystem_2D' : 2, 'shallow_roe_with_efix_2D' : 3, + 'shallow_bathymetry_fwave_2D' : 3, + 'sw_aug_2D' : 3, 'shallow_sphere_2D' : 3, 'vc_advection_2D' : 1, 'vc_acoustics_2D' : 2,