diff --git a/examples/bouss/README.txt b/examples/bouss/README.txt index d490d3591..32e0b6f80 100644 --- a/examples/bouss/README.txt +++ b/examples/bouss/README.txt @@ -28,9 +28,6 @@ variables properly, see: The file setenv.sh illustrates how you might set some environment variables for the bash shell. -The file setenv.sh illustrates how you might set some environment -variables for the bash shell. - A file petscMPIoptions is also required to set some PETSc parameters for the iterative method used to solve the large sparse linear systems that arise at each refinement level when the Boussinesq equations are solved. diff --git a/examples/storm-surge/isaac/test_isaac.py b/examples/storm-surge/isaac/test_isaac.py index 9c5585df4..65adbfb0e 100644 --- a/examples/storm-surge/isaac/test_isaac.py +++ b/examples/storm-surge/isaac/test_isaac.py @@ -4,16 +4,15 @@ Primarily tests all of the input formats GeoClaw handles """ -from pathlib import Path -import os import gzip +import os import unittest -import pytest - -import numpy as np +from pathlib import Path -import clawpack.geoclaw.test as test import clawpack.clawutil as clawutil +import clawpack.geoclaw.test as test +import numpy as np +import pytest from clawpack.geoclaw.surge.storm import Storm days2seconds = lambda t: t * 60.0**2 * 24.0 diff --git a/src/2d/bouss/Makefile.bouss b/src/2d/bouss/Makefile.bouss index 0fa48b96f..256b102c3 100644 --- a/src/2d/bouss/Makefile.bouss +++ b/src/2d/bouss/Makefile.bouss @@ -159,6 +159,7 @@ COMMON_SOURCES += \ $(BOUSSLIB)/buildSparseMatrixSGNcrs.f90 \ $(BOUSSLIB)/prepBuildSparseMatrixSGNcoo.f90 \ $(BOUSSLIB)/buildSparseMatrixSGNcoo.f90 \ + $(BOUSSLIB)/buildSparseMatrixSGNcoo_blocks.f90 \ $(BOUSSLIB)/buildSparseMatrixMScoo.f90 \ $(BOUSSLIB)/compressOut.f \ $(BOUSSLIB)/matvec.f90 \ diff --git a/src/2d/bouss/amr2.f90 b/src/2d/bouss/amr2.f90 index 09e2b241b..6cc9bacde 100644 --- a/src/2d/bouss/amr2.f90 +++ b/src/2d/bouss/amr2.f90 @@ -860,7 +860,6 @@ program amr2 write(*,format_string) write(*,*) write(timing_unit,*) - close(timing_unit) if (isolver.eq. 2) then ! if (countIterRef == 0) then @@ -933,6 +932,7 @@ program amr2 close(dbugunit) #ifdef HAVE_PETSC + call PetscViewerASCIIStdoutSetFileUnit(timing_unit,ierr) call PetscFinalize(ierr) #endif diff --git a/src/2d/bouss/bouss_module.f90 b/src/2d/bouss/bouss_module.f90 index 4cf8cc705..8ff779f92 100644 --- a/src/2d/bouss/bouss_module.f90 +++ b/src/2d/bouss/bouss_module.f90 @@ -19,7 +19,7 @@ module bouss_module logical :: startWithBouss - logical :: crs + logical :: crs, origCooFormat real(kind=8) :: startBoussTime @@ -165,8 +165,17 @@ subroutine set_bouss(rest,time,naux) ! using compressed row storage (CRS, also known as CSR or Yale) format: crs = .true. endif - - ! crs = .false. ! uncomment this line to force COO with SGN for testing + ! if crs is false this will force SGN to use triplet COO form + ! crs = .false. + + ! origCoo format ordered unknowns by all u updates then all v updates + ! block format is u and v together for a given cell. Better cache + ! performance, better for debugging and comparing with crs. + ! also reads into matlab for cond number testing, plotting, + ! so keeping around + !origCooFormat = .false. ! set to false to use block format with coo + ! setting to avoid uninitialized var, but not used unless crs set to false + origCooFormat = .true. !------------------------------------------ if (rest) then diff --git a/src/2d/bouss/buildSparseMatrixSGNcoo.f90 b/src/2d/bouss/buildSparseMatrixSGNcoo.f90 index 4649f9266..fcf71a1eb 100644 --- a/src/2d/bouss/buildSparseMatrixSGNcoo.f90 +++ b/src/2d/bouss/buildSparseMatrixSGNcoo.f90 @@ -952,7 +952,7 @@ subroutine buildSparseMatrixSGNcoo(q,qold,aux,soln,rhs, & ! Need to add comments on SGN RHS. if (debug) then - write(67,*) '++b i,j,rhs = ',i,j,rhs(k_ij),rhs(k_ij + nD) + write(67,*) '++before i,j,rhs = ',i,j,rhs(k_ij),rhs(k_ij + nD) endif !RHS for SGN @@ -976,20 +976,26 @@ subroutine buildSparseMatrixSGNcoo(q,qold,aux,soln,rhs, & + 0.5d0*h_ij*wy + w(i,j)*etay) + if (debug) then + write(67,*) '++after i,j,rhs = ',i,j,rhs(k_ij),rhs(k_ij + nD) + endif + if (debug) then write(44,998) nelt,nelt-neltSt 998 format(" ending with nelt ",i6," used ",i5) endif - enddo - enddo - if (debug .and. .false.) then ! dump matrix to look for singularity - write(88,*)" level ",levelBouss + if (debug) then ! dump matrix to look for singularity + write(88,*)" triplets for i,j ",i,j," level ",levelBouss do k = neltBegin+1, nelt - write(88,103) minfo_matrix_ia(k),minfo_matrix_ja(k),minfo_matrix_sa(k) - 103 format(2i8,e16.7) + !write(88,103) minfo_matrix_ia(k),minfo_matrix_ja(k),minfo_matrix_sa(k) + write(88,103) minfo_matrix_sa(k) + !103 format(2i8,e16.7) + 103 format(e16.7) end do + write(88,*) + !close(88) !write(89,*)" level ",levelBouss !do k = 1,2*numBoussCells @@ -1000,6 +1006,9 @@ subroutine buildSparseMatrixSGNcoo(q,qold,aux,soln,rhs, & !!stop endif + enddo + enddo + if (debug) then if (mptr .eq. 1) then write(88,*)" grid ",mptr," level ",levelBouss,"x,y,ia,ja,sa" diff --git a/src/2d/bouss/buildSparseMatrixSGNcoo_blocks.f90 b/src/2d/bouss/buildSparseMatrixSGNcoo_blocks.f90 new file mode 100644 index 000000000..6a67e0bff --- /dev/null +++ b/src/2d/bouss/buildSparseMatrixSGNcoo_blocks.f90 @@ -0,0 +1,1117 @@ + +subroutine buildSparseMatrixSGNcoo_blocks(q,qold,aux,soln,rhs, & + minfo_matrix_ia,minfo_matrix_ja,minfo_matrix_sa, & + numBoussCells,levelBouss, & + mptr,nx,ny,nvar,naux, & + matrix_rhs1_max,matrix_rhs2_max) + + + use amr_module + use geoclaw_module, only: earth_radius, deg2rad, coordinate_system, grav + use geoclaw_module, only: sea_level, dry_tolerance + use bouss_module + + implicit none + + integer, intent(in) :: nvar, naux, levelBouss, numBoussCells + + !! pass in numBoussCells for this level so can dimension these arrays + real(kind=8) :: soln(0:2*numBoussCells), rhs(0:2*numBoussCells) + + type(matrix_patchIndex), pointer :: mi + type(matrix_levInfo), pointer :: minfo + + integer :: i,j,nng,nb,nelt,levSt,nx,ny,loc,locaux,mptr + integer :: mitot, mjtot, neltSt, neltBegin + integer :: jst, jend, ivar + + real(kind=8) :: q(nvar,1-nghost:nx+nghost,1-nghost:ny+nghost) + real(kind=8) :: qold(nvar,1-nghost:nx+nghost,1-nghost:ny+nghost) + real(kind=8) :: aux(naux,1-nghost:nx+nghost,1-nghost:ny+nghost) + + real(kind=8), dimension(0:nx+1,0:ny+1) :: Bxx, Byy, Bx, By, Bxy + real(kind=8), dimension(0:nx+1,0:ny+1) :: phi, eta, w + real(kind=8) :: etax, etay, hx,hy, wx,wy + real(kind=8) :: xlow, ylow, x, y + integer :: ii,jj,kd + + integer :: nD + integer :: k, k_ij, k_imj, k_ipj, k_ijm, k_ijp + integer :: m,n + integer :: k_imjp, k_imjm, k_ipjm, k_ipjp + integer :: kh_ij, kh_imj, kh_ipj, kh_ijm, kh_ijp + integer :: kh_imjp, kh_imjm, kh_ipjm, kh_ipjp, khu_ij, khv_ij + + real(kind=8) :: h_ij + real(kind=8) :: dx, dy, dxm, dym, dxdx, dydy, dxdy4 + + integer :: nelt_ij, nelt_ij_nD + real(kind=8) :: sa, huc, hvc, matrix_rhs1_max, matrix_rhs2_max + integer :: ja, bc, c + real(kind=8) :: b_min, b_max + integer :: minfo_matrix_ia(24*numBoussCells) + integer :: minfo_matrix_ja(24*numBoussCells) + real(kind=8) :: minfo_matrix_sa(24*numBoussCells) + real(kind=8) :: u(1-nghost:nx+nghost,1-nghost:ny+nghost), v(1-nghost:nx+nghost,1-nghost:ny+nghost) + real(kind=8) :: phix, phiy + + logical :: debug + logical :: yhi_db, ylo_db, xlo_db, xhi_db + logical :: allzero + + + ! bc_xlo, bc_ylo, bc_xhi, bc_yhi are set in bouss_module. + ! 0 means Dirchlet value 0 (no correction at each of union of patches) + ! 1 means Neumann (correction in first interior cell = in ghost cell) + ! 2 means use ghost cell values in q(4:5,:,:) as Dirichlet BCs + ! Multiply by matrix element and move to RHS + ! Requires nvar==5. + + +#ifdef WHERE_AM_I + write(*,*) "starting buildSparseMatrixSGNcoo" +#endif + + debug = .false. + !debug = .true. + !write(67,*)"level ",levelBouss + + minfo => matrix_info_allLevs(levelBouss) + + !write(31,*)"Start of buildSparseMatrixSGNcoo for grid ",mptr + !write(*,*)"Start of buildSparseMatrixSGNcoo for grid ",mptr + mitot = nx + 2*nghost + mjtot = ny + 2*nghost + xlow = rnode(cornxlo,mptr) + ylow = rnode(cornylo,mptr) + !call symcheck2(q,mitot,mjtot) + + +! +! ================== Step 2 set up for matrix elements ====================== + + ! Set up matrix A and RHS by looping over all Bouss patches: + + if (debug) then + write(49,*)"Level ",levelBouss," grid ",mptr," at start of buildSparseMatrixSGNcoo" + do j = 1-nghost, ny+nghost + do i = 1-nghost, nx+nghost + write(49,900)i+nghost,j+nghost,(q(m,i,j),m=1,nvar) + 900 format(2i5,5e15.6) + end do + end do + endif + + dx = hxposs(levelBouss) + dy = hyposs(levelBouss) + + + if (coordinate_system .eq. 2) then + ! convert from dy in degrees to meters: + dym = dy * earth_radius * deg2rad + ! dxm depends on latitude, done in loops below instead of here + else + dym = dy + endif + + dydy = dym**2 + nD = numBoussCells ! shorthand for size of matrix blocks D1 to D4 + nelt = minfo%matrix_nelt ! some number of matrix elements already added, + neltBegin = nelt ! this is where this grid started + + !REORG to put u and v together in blocks + + ! set booleans indicating whether each patch edge is a Domain Boundary: + !call setDomainBndries(mptr,yhi_db,ylo_db,xlo_db,xhi_db) + xlo_db = (node(ndilo,mptr) .eq. 0) + xhi_db = (node(ndihi,mptr) .eq. (iregsz(levelBouss)-1)) + ylo_db = (node(ndjlo,mptr) .eq. 0) + yhi_db = (node(ndjhi,mptr) .eq. (jregsz(levelBouss)-1)) + + + + ! Get integer array providing matrix indices mapping + ! cell coordinates i,j to location k in vector of all unknowns: + + nb = node(bpatchNum,mptr) + mi => minfo%matrix_indices(nb) ! integer array for this patch + + ! temp arrays needed for this patch: + ! provide one layer of ghost cells needed when computing + ! second differences of these quantities, but values in ghost + ! cells are never actually needed since second diff's won't be + ! needed in first row of interior cells + + ! First loop over all cells (i,j), one row of ghost cells, + ! computing s1, s2, h0etax, h0etay, + ! quantities that involve first derivatives. In the second loop + ! over all cells below we will compute matrix_rhs using second + ! differences of these quantities. + + ! compute eta used below in update + ! Remember that ETA not same size as Q + eta(:,:) = q(1,0:nx+1,0:ny+1) + & + aux(1,0:nx+1,0:ny+1) + + ! will have to FIX for h = 0? Hopefully not that close to 0 + !u(:,:) = q(2,:,:) / q(1,:,:) + !v(:,:) = q(3,:,:) / q(1,:,:) + !where (q(1,:,:) /= 0.d0) + where (abs(q(1,:,:)) > dry_tolerance) + u(:,:) = q(2,:,:) / q(1,:,:) + v(:,:) = q(3,:,:) / q(1,:,:) + elsewhere + u(:,:) = 0.d0 + v(:,:) = 0.d0 + end where + + + do j=0,ny+1 + if (coordinate_system .eq. 2) then + ! aux(3,i,j) is latitude, doesn't depend on i so use i=1 + dxm = dx * earth_radius * deg2rad * cos(aux(3,1,j)) + else + dxm = dx + endif + dxdx = dxm**2 + dxdy4 = 4.d0*dxm*dym + + ! compute phi. rest will now be computed ! in second loop + do i=0,nx+1 + + phi(i,j) = ((v(i+1,j)-v(i-1,j))*(u(i,j+1)-u(i,j-1)) - & + (u(i+1,j)-u(i-1,j))*(v(i,j+1)-v(i,j-1)))/dxdy4 & + + ((u(i+1,j)-u(i-1,j))/(2.d0*dxm) & + + (v(i,j+1)-v(i,j-1))/(2.d0*dym))**2 + Bx(i,j) = (aux(1,i+1,j) - aux(1,i-1,j))/(2.d0*dxm) + By(i,j) = (aux(1,i,j+1) - aux(1,i,j-1))/(2.d0*dym) + Bxy(i,j) = (aux(1,i+1,j+1)-aux(1,i+1,j-1) & + -aux(1,i-1,j+1)+aux(1,i-1,j-1))/dxdy4 + Bxx(i,j) = (aux(1,i+1,j)-2.d0*aux(1,i,j)+aux(1,i-1,j))/dxm**2 + Byy(i,j) = (aux(1,i,j+1)-2.d0*aux(1,i,j)+aux(1,i,j-1))/dym**2 + w(i,j) = u(i,j)**2 * Bxx(i,j) + 2.d0*u(i,j)*v(i,j)*Bxy(i,j) & + +v(i,j)**2 * Byy(i,j) + + enddo + enddo + +! ================== build matrix and rhs + ! Second loop over all cells. ONLY INTERIOR OK NOW + + do j=1,ny + if (coordinate_system .eq. 2) then + ! aux(3,i,j) is latitude, doesn't depend on i so use i=1 + dxm = dx * earth_radius * deg2rad * & + cos(aux(3,1,j)) + else + dxm = dx + endif + dxdx = dxm**2 + dxdy4 = 4.d0*dxm*dym + + + do i=1,nx + + ! sparse matrix row indices for cell (i,j) are: + ! k_ij for the row of I-D11 and -D12 in block matrix, + ! k_ij + nD for the row of -D21 and I-D22 + + ! in this inner loop we set all nonzero matrix elements + ! in these two rows of the matrix (and rhs) corresponding + ! to momenumtum updates to hu and hv in the (i,j) cell + + ! sparse matrix column indices for neighboring cells + ! are the values k computed next (for I-D1 and -D3) + ! and the corresponding k + nD for the corresponding column + ! of the -D2 and I-D4. These indices are needed to + ! set matrix elements in these columns, corresponding to + ! momentum components hu and hv in the neighboring cells + + k_imj = mi%mindex(i-1,j) + k_ipj = mi%mindex(i+1,j) + k_ijm = mi%mindex(i,j-1) + k_ijp = mi%mindex(i,j+1) + k_imjm = mi%mindex(i-1,j-1) + k_imjp = mi%mindex(i-1,j+1) + k_ipjm = mi%mindex(i+1,j-1) + k_ipjp = mi%mindex(i+1,j+1) + + ! use the fact that h,hu,hv are contiguous in alloc: + h_ij = q(1,i,j) ! to save typing keep this one + + k_ij = mi%mindex(i,j) + if (k_ij == -1) then + write(6,*) '+++ unexpected k_ij=-1 at i,j: ',i,j + endif + neltSt = nelt + if (debug) then + write(44,999)i,j,k_ij,k_ij,nelt + 999 format("cell ",2i4," k_ij = ",i5," before shift was ",i5," starting with nelt ",i6) + endif + + !b_min = min(aux(1,i,j), aux(1,i-1,j),aux(1,i+1,j), & + ! aux(1,i,j-1),aux(1,i,j+1), & + ! aux(1,i-1,j-1),aux(1,i-1,j+1), & + ! aux(1,i+1,j-1),aux(1,i+1,j+1)) + !b_max = max(aux(1,i,j), aux(1,i-1,j),aux(1,i+1,j), & + ! aux(1,i,j-1),aux(1,i,j+1), & + ! aux(1,i-1,j-1),aux(1,i-1,j+1), & + ! aux(1,i+1,j-1),aux(1,i+1,j+1)) + if (.not. mi%isBouss(i,j)) then + + ! The deepest water in block of 9 cells is shallower + ! than boussMinDepth, or else one of the cells is onshore. + + ! Revert to SWE in this cell, i.e. two rows + ! of matrix should be rows of identity with rhs = 0, + ! so correction terms huc,hvc calculated will be 0. + + !QUESTION: I should be h_ij*I. FIXED + ! BUT: then discovered the leading h in matrix is a typo + ! fixed for now by multiplying rhs by h_ij too. + + ! Replace diagonal block I-D11 by I matrix: + nelt = nelt + 1 + !minfo_matrix_ia(nelt) = k_ij + !minfo_matrix_ja(nelt) = k_ij + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = 2*k_ij-1 + + ! changed to +1*max on sept 2 2022 to test gamg + minfo_matrix_sa(nelt) = 1.d0*max(abs(h_ij),1000.d0) + + ! Replace diagonal block I-D22 by I matrix: + nelt = nelt + 1 + !minfo_matrix_ia(nelt) = k_ij + nD + !minfo_matrix_ja(nelt) = k_ij + nD + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = 2*k_ij + + ! changed to +1*max on sept 2 2022 to test gamg + minfo_matrix_sa(nelt) = 1.d0*max(abs(h_ij),1000.d0) + + ! -D12 and -D21 are the 0 matrix in this case (no nozeros) + + cycle ! go on to next i + endif + + + ! If we get here, we are dealing with interior point + ! Add rows k_ij and (k_ij + nD) to matrix: + ! Row k_ij corresponds to correction to hu in cell (i,j) + ! Row k_ij+nD corresponds to correction to hv in cell (i,j) + ! The nonzero matrix elements will be in the same columns + ! and also in columns + ! k_ipj, k_ipj+nD for the coeff of the (i+1,j) stencil pt + ! k_imjp, k_imjp+nD for the coeff of (i-1,j+1) point, + ! etc. + ! Note that if e.g. k_ipj = -1 then cell (i+1,j) is not in the + ! union of patches at this level and so there is no + ! corresponding column. + ! In this case: + ! For Dirichlet BC with value 0, the stencil value is discarded, + ! For Neumann BC, the stencil coefficient is instead applied + ! to the neighboring cell that is in the patch, + ! assuming constant extrapolation at each edge. + + ! Increment nelt each time we add another nonzero matrix + ! element. + ! Store index nelt_ij for index to diagonal element in row k_ij + ! Store index nelt_ij_nD for diagonal element in row k_ij+nD + + ! initialize to -1 for debugging in case not set properly: + nelt_ij = -1 + nelt_ij_nD = -1 + + hx = (q(1,i+1,j) - q(1,i-1,j))/(2.d0*dxm) + hy = (q(1,i,j+1) - q(1,i,j-1))/(2.d0*dym) + etax = (eta(i+1,j) - eta(i-1,j))/(2.d0*dxm) + etay = (eta(i,j+1) - eta(i,j-1))/(2.d0*dym) + + !----------------------------------------------- + ! I-D11 block (x and xx derivatives of hu components): + + ! Matrix element row U(i,j) column U(i,j): + sa = (1.d0 + alpha*(2.d0*h_ij**2/(3.d0*dxdx) & + + 0.5d0*h_ij*Bxx(i,j) + Bx(i,j)*etax)) + nelt = nelt+1 + !minfo_matrix_ia(nelt) = k_ij + !minfo_matrix_ja(nelt) = k_ij + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = 2*k_ij-1 + minfo_matrix_sa(nelt) = sa + + nelt_ij = nelt ! index of diagonal matrix element + + + ! Matrix element row U(i,j) column U(i+1,j): + sa = alpha*(-h_ij*hx/(2.d0*dxm) - h_ij**2/(3.d0*dxdx)) + + if (k_ipj > -1) then + ! adjacent cell is interior pt, set matrix element to sa: + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = 2*k_ipj-1 + minfo_matrix_sa(nelt) = sa + else + ! right boundary of patch + if (.not. xhi_db) then + ! Dirichlet BC using Bouss correction from ghost cell + huc = q(4,i+1,j) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*huc + else if (bc_xhi==1) then + ! For Neumann BC add sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij) = minfo_matrix_sa(nelt_ij) + sa + else if (bc_xhi==3) then + ! For wall, negate u: add -sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij) = minfo_matrix_sa(nelt_ij) - sa + endif + endif + + ! Matrix element row U(i,j) column U(i-1,j): + sa = alpha*(h_ij*hx/(2.d0*dxm) - h_ij**2/(3.d0*dxdx)) + if (k_imj > -1) then + ! adjacent cell is interior pt, set matrix element to sa: + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = 2*k_imj-1 + minfo_matrix_sa(nelt) = sa + else + ! left boundary of patch + if (.not. xlo_db) then + ! Dirichlet BC using Bouss correction from ghost cell + huc = q(4,i-1,j) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*huc + else if (bc_xlo==1) then + ! For Neumann BC add sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij) = minfo_matrix_sa(nelt_ij) + sa + else if (bc_xlo==3) then + ! For wall, negate u: add -sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij) = minfo_matrix_sa(nelt_ij) - sa + endif + endif + + + !----------------------------------------------- + ! I-D22 block (yy derivatives of hv components): + + ! Matrix element row V(i,j) column V(i,j): + sa = (1.d0 + alpha*(2.d0*h_ij**2/(3.d0*dydy) & + + 0.5d0*h_ij*Byy(i,j) + By(i,j)*etay)) + nelt = nelt+1 + !minfo_matrix_ia(nelt) = k_ij + nD + !minfo_matrix_ja(nelt) = k_ij + nD + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = 2*k_ij + minfo_matrix_sa(nelt) = sa + + nelt_ij_nD = nelt ! index of diagonal matrix element + + + ! Matrix element row V(i,j) column V(i,j+1): + sa = alpha*(-h_ij*hy/(2.d0*dym) - h_ij**2/(3.d0*dydy)) + if (k_ijp > -1) then + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = 2*k_ijp + minfo_matrix_sa(nelt) = sa + else + ! top boundary of patch + if (.not. yhi_db) then + ! Dirichlet BC using Bouss correction from ghost cell + hvc = q(5,i,j+1) + !rhs(k_ij+nD) = rhs(k_ij+nD) - sa*hvc + rhs(2*k_ij) = rhs(2*k_ij) - sa*hvc + else if (bc_yhi==1) then + ! For Neumann BC add sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij_nD) = & + minfo_matrix_sa(nelt_ij_nD) + sa + else if (bc_yhi==3) then + ! For wall BC add -sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij_nD) = & + minfo_matrix_sa(nelt_ij_nD) - sa + endif + endif + + + ! Matrix element row V(i,j) column V(i,j-1): + sa = alpha*(h_ij*hy/(2.d0*dym) - h_ij**2/(3.d0*dydy)) + if (k_ijm > -1) then + nelt = nelt+1 + !minfo_matrix_ia(nelt) = k_ij + nD + !minfo_matrix_ja(nelt) = k_ijm + nD + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = 2*k_ijm + minfo_matrix_sa(nelt) = sa + else + ! bottom boundary of patch + if (.not. ylo_db) then + ! Dirichlet BC using Bouss correction from ghost cell + hvc = q(5,i,j-1) + !rhs(k_ij+nD) = rhs(k_ij+nD) - sa*hvc + rhs(2*k_ij) = rhs(2*k_ij) - sa*hvc + else if (bc_ylo==1) then + ! For Neumann BC add sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij_nD) = & + minfo_matrix_sa(nelt_ij_nD) + sa + else if (bc_ylo==3) then + ! For wall BC add -sa to diagonal matrix element: + minfo_matrix_sa(nelt_ij_nD) = & + minfo_matrix_sa(nelt_ij_nD) - sa + endif + endif + + + !----------------------------------------------- + ! -D12 block (cross derivatives): + ! row ia = k_ij for hu corrections in cell (i,j) + ! columns k_ipjp + nD and other neighbors for hv corrections + + ! Matrix element row U(i,j) column V(i,j): + ! NEW MATRIX ELEMENT FOR SGN WITH VARYING TOPO + sa = alpha*(0.5d0*h_ij*Bxy(i,j) + etax*By(i,j)) + nelt = nelt+1 + !minfo_matrix_ia(nelt) = k_ij + !minfo_matrix_ja(nelt) = k_ij + nD + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = 2*k_ij + minfo_matrix_sa(nelt) = sa + + + ! Matrix element row U(i,j) column V(i,j+1): + ! NEW MATRIX ELEMENT FOR SGN + sa = -alpha*h_ij*(hx + 0.5d0*Bx(i,j))/(2.d0*dym) + ja = -1 + if (k_ijp > -1) then + !ja = k_ijp + nD + ja = 2*k_ijp + else if (.not. yhi_db) then + hvc = q(5,i,j+1) + !rhs(k_ij) = rhs(k_ij) - sa*hvc + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else + ! at top boundary, use cell ij instead, + !ja = k_ij + nD + ja = 2*k_ij + if (bc_yhi==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i,j-1): + ! NEW MATRIX ELEMENT FOR SGN + sa = alpha*(h_ij*(hx + 0.5d0*Bx(i,j))/(2.d0*dym)) + ja = -1 + if (k_ijm > -1) then + !ja = k_ijm + nD + ja = 2*k_ijm + else if (.not. ylo_db) then + hvc = q(5,i,j-1) + !rhs(k_ij) = rhs(k_ij) - sa*hvc + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else + ! at bottom boundary, use cell ij instead, + !ja = k_ij + nD + ja = 2*k_ij + if (bc_ylo==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i+1,j): + ! NEW MATRIX ELEMENT FOR SGN WITH TOPO + sa = alpha*0.5d0*h_ij*By(i,j) / (2.d0*dxm) + ja = -1 + if (k_ipj > -1) then + !ja = k_ipj + nD + ja = 2*k_ipj + else if (.not. xhi_db) then + hvc = q(5,i+1,j) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else + ! at right boundary, use cell ij instead, + !ja = k_ij + nD + ja = 2*k_ij + ! no reflection of V at right wall + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i-1,j): + ! NEW MATRIX ELEMENT FOR SGN WITH TOPO + sa = -alpha*0.5d0*h_ij*By(i,j) / (2.d0*dxm) + ja = -1 + if (k_imj > -1) then + !ja = k_imj + nD + ja = 2*k_imj + else if (.not. xlo_db) then + hvc = q(5,i-1,j) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else + ! at left boundary, use cell ij instead, + !ja = k_ij + nD + ja = 2*k_ij + ! no reflection of V at left wall + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i+1,j+1): + sa = -alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_ipjp > -1) then + !ja = k_ipjp + nD + ja = 2*k_ipjp + !else if (((k_ijp == -1) .and. (.not. yhi_db)) .or. & + ! ((k_ipj == -1) .and. (.not. xhi_db))) then + else if ((.not. yhi_db .or. j .lt. ny) .and. & + (.not. xhi_db .or. i .lt. nx)) then + hvc = q(5,i+1,j+1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else if (k_ijp > -1) then + ! at right boundary, use cell above instead + !ja = k_ijp + nD + ja = 2*k_ijp + else if (k_ipj > -1) then + ! at top boundary, use cell to right instead, + !ja = k_ipj + nD + ja = 2*k_ipj + if (bc_yhi==3) sa = -sa ! reflect V in y + else if (xhi_db .and. i .eq. nx .and. .not. yhi_db) then + !ja = k_ij + nD + ja = 2*k_ij + else + ! top-right corner + !ja = k_ij + nD + ja = 2*k_ij + if (bc_yhi==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i+1,j-1): + sa = alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_ipjm > -1) then + !ja = k_ipjm + nD + ja = 2*k_ipjm + !else if (((k_ijm == -1) .and. (.not. ylo_db)) .or. & + ! ((k_ipj == -1) .and. (.not. xhi_db))) then + else if ((.not. ylo_db .or. j>1) .and. & + (.not. xhi_db .or. i .lt. nx)) then + hvc = q(5,i+1,j-1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else if (k_ijm > -1) then + ! at right boundary, use cell below instead + !ja = k_ijm + nD + ja = 2*k_ijm + else if (k_ipj > -1) then + ! at bottom boundary, use cell to right instead, + !ja = k_ipj + nD + ja = 2*k_ipj + if (bc_ylo==3) sa = -sa ! reflect V in y + else if (xhi_db .and. i.eq.nx .and. .not. ylo_db) then + !ja = k_ij + nD + ja = 2*k_ij + else + ! bottom-right corner + !ja = k_ij + nD + ja = 2*k_ij + if (bc_ylo==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row U(i,j) column V(i-1,j+1): + sa = alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_imjp > -1) then + !ja = k_imjp + nD + ja = 2*k_imjp + !else if (((k_ijp == -1) .and. (.not. yhi_db)) .or. & + ! ((k_imj == -1) .and. (.not. xlo_db))) then + else if ((.not. yhi_db .or. j .lt. ny) .and. & + (.not. xlo_db .or. i > 1)) then + hvc = q(5,i-1,j+1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else if (k_ijp > -1) then + ! at left domain boundary, use cell above instead + !ja = k_ijp + nD + ja = 2*k_ijp + else if (k_imj > -1) then + ! at top domain boundary, use cell to left instead, + !ja = k_imj + nD + ja = 2*k_imj + if (bc_yhi==3) sa = -sa ! reflect V in y + else if (xlo_db .and. i .eq. 1 .and. .not. yhi_db) then + !ja = k_ij + nD + ja = 2*k_ij + else + ! top-left corner + !ja = k_ij + nD + ja = 2*k_ij + if (bc_yhi==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + ! Matrix element row U(i,j) column V(i-1,j-1): + sa = -alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_imjm > -1) then + ja = 2*k_imjm + !else if (((k_ijm == -1) .and. (.not. ylo_db)) .or. & + ! ((k_imj == -1) .and. (.not. xlo_db))) then + else if ((.not. ylo_db .or. j>1) .and. & + (.not. xlo_db .or. i>1)) then + hvc = q(5,i-1,j-1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*hvc + else if (k_ijm > -1) then + ! at left boundary, use cell below instead + !ja = k_ijm + nD + ja = 2*k_ijm + else if (k_imj > -1) then + ! at bottom boundary, use cell to left instead, + !ja = k_imj + nD + ja = 2*k_imj + if (bc_ylo==3) sa = -sa ! reflect V in y + else if (xlo_db .and. i .eq. 1 .and. .not. ylo_db) then + !ja = k_ij + nD + ja = 2*k_ij + else + ! bottom-left corner + !ja = k_ij + nD + ja = 2*k_ij + if (bc_ylo==3) sa = -sa ! reflect V in y + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij-1 + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + !----------------------------------------------- + ! -D21 block (cross derivatives + ! row ia = k_ij + nD for hv corrections in cell (i,j) + ! columns k_ipjp and other neighbors for hu corrections + + ! Matrix element row V(i,j) column U(i,j): + ! NEW MATRIX ELEMENT FOR SGN WITH VARYING TOPO + sa = alpha*(0.5d0*h_ij*Bxy(i,j) + etay*Bx(i,j)) + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = 2*k_ij-1 + minfo_matrix_sa(nelt) = sa + + ! Matrix element row V(i,j) column U(i+1,j): + ! NEW MATRIX ELEMENT FOR SGN + sa = -alpha*h_ij*(hy + 0.5d0*By(i,j))/(2.d0*dxm) + ja = -1 + if (k_ipj > -1) then + !ja = k_ipj + ja = 2*k_ipj-1 + else if (.not. xhi_db) then + huc = q(4,i+1,j) + !rhs(k_ij+nD) = rhs(k_ij+nD) - sa*huc + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else + ! at right boundary, use cell ij instead, + !ja = k_ij + ja = 2*k_ij-1 + if (bc_xhi==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i-1,j): + ! NEW MATRIX ELEMENT FOR SGN + sa = alpha*h_ij*(hy + 0.5d0*By(i,j))/(2.d0*dxm) + ja = -1 + if (k_imj > -1) then + !ja = k_imj + ja = 2*k_imj-1 + else if (.not. xlo_db) then + huc = q(4,i-1,j) + !rhs(k_ij+nD) = rhs(k_ij+nD) - sa*huc + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else + ! at left boundary, use cell ij instead, + !ja = k_ij + ja = 2*k_ij-1 + if (bc_xlo==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i,j+1): + ! NEW MATRIX ELEMENT FOR SGN WITH TOPO + sa = alpha*0.5d0*h_ij*Bx(i,j) / (2.d0*dym) + ja = -1 + if (k_ijp > -1) then + !ja = k_ijp ! FIXED BUG2 + ja = 2*k_ijp-1 ! FIXED BUG2 + else if (.not. yhi_db) then + huc = q(4,i,j+1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*huc + else + ! at top boundary, use cell ij instead, + !ja = k_ij + ja = 2*k_ij-1 + ! no reflection of U at top wall + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i,j-1): ! FIXED BUG2 + ! NEW MATRIX ELEMENT FOR SGN WITH TOPO + sa = -alpha*0.5d0*h_ij*Bx(i,j) / (2.d0*dym) + ja = -1 + if (k_ijm > -1) then + !ja = k_ijm ! FIXED BUG2 + ja = 2*k_ijm-1 ! FIXED BUG2 + else if (.not. ylo_db) then + huc = q(4,i,j-1) + rhs(2*k_ij-1) = rhs(2*k_ij-1) - sa*huc + else + ! at bottom boundary, use cell ij instead, + !ja = k_ij ! FIXED BUG2 + ja = 2*k_ij-1 ! FIXED BUG2 + ! no reflection of U at bottom wall + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i+1,j+1): + sa = -alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_ipjp > -1) then + !ja = k_ipjp + ja = 2*k_ipjp-1 + !else if (((k_ijp == -1) .and. (.not. yhi_db)) .or. & + ! ((k_ipj == -1) .and. (.not. xhi_db))) then + else if ((.not. yhi_db .or. j.lt.ny) .and. & + (.not. xhi_db .or. i.lt.nx)) then + huc = q(4,i+1,j+1) + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else if (k_ijp > -1) then + ! at right boundary, use cell above instead, + !ja = k_ijp + ja = 2*k_ijp-1 + if (bc_xhi==3) sa = -sa ! reflect U in x + else if (k_ipj > -1) then + ! at top boundary, use cell to right instead, + !ja = k_ipj + ja = 2*k_ipj-1 + else if (yhi_db .and. j .eq. ny .and. .not. xhi_db) then + !ja = k_ij + ja = 2*k_ij-1 + else + ! top-right corner + !ja = k_ij ! FIXED BUG + ja = 2*k_ij-1 ! FIXED BUG + if (bc_xhi==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i+1,j-1): + sa = alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_ipjm > -1) then + !ja = k_ipjm + ja = 2*k_ipjm-1 + !else if (((k_ijm == -1) .and. (.not. ylo_db)) .or. & + ! ((k_ipj == -1) .and. (.not. xhi_db))) then + else if ((.not. ylo_db .or. j>1) .and. & + (.not. xhi_db .or. i .lt. nx)) then + huc = q(4,i+1,j-1) + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else if (k_ijm > -1) then + ! at right boundary, use cell below instead, + !ja = k_ijm + ja = 2*k_ijm-1 + if (bc_xhi==3) sa = -sa ! reflect U in x + else if (k_ipj > -1) then + ! at bottom boundary, use cell to right instead, + !ja = k_ipj + ja = 2*k_ipj-1 + else if (ylo_db .and. j .eq. 1 .and. .not. xhi_db) then + !ja = k_ij + ja = 2*k_ij-1 + else + ! bottom-right corner + !ja = k_ij ! FIXED BUG + ja = 2*k_ij-1 ! FIXED BUG + if (bc_xhi==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + if (debug) write(67,*) i,j,i+1,j-1,c,nelt,sa + endif + + + ! Matrix element row V(i,j) column U(i-1,j+1): + sa = alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_imjp > -1) then + !ja = k_imjp + ja = 2*k_imjp-1 + !else if (((k_ijp == -1) .and. (.not. yhi_db)) .or. & + ! ((k_imj == -1) .and. (.not. xlo_db))) then + else if ((.not. yhi_db .or. j .lt. ny) .and. & + (.not. xlo_db .or. i>1)) then + huc = q(4,i-1,j+1) + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else if (k_ijp > -1) then + ! at left boundary, use cell above instead, + !ja = k_ijp + ja = 2*k_ijp-1 + if (bc_xlo==3) sa = -sa ! reflect U in x + else if (k_imj > -1) then + ! at top boundary, use cell to left instead, + !ja = k_imj + ja = 2*k_imj-1 + else if (yhi_db .and. j.eq.ny .and. .not. xlo_db) then + !ja = k_ij + ja = 2*k_ij-1 + else + ! top-left corner + !ja = k_ij !FIXED BUG + ja = 2*k_ij-1 !FIXED BUG + if (bc_xlo==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + ! Matrix element row V(i,j) column U(i-1,j-1): + sa = -alpha*h_ij**2 / (3.d0*dxdy4) + ja = -1 + if (k_imjm > -1) then + !ja = k_imjm + ja = 2*k_imjm-1 + !else if (((k_ijm == -1) .and. (.not. ylo_db)) .or. & + ! ((k_imj == -1) .and. (.not. xlo_db))) then + else if ((.not. ylo_db .or. j>1) .and. & + (.not. xlo_db .or. i>1)) then + huc = q(4,i-1,j-1) + rhs(2*k_ij) = rhs(2*k_ij) - sa*huc + else if (k_ijm > -1) then + ! at left boundary, use cell below instead, + !ja = k_ijm + ja = 2*k_ijm-1 + if (bc_xlo==3) sa = -sa ! reflect U in x + else if (k_imj > -1) then + ! at bottom boundary, use cell to left instead, + !ja = k_imj + ja = 2*k_imj-1 + else if (ylo_db .and. j .eq. 1 .and. .not. xlo_db) then + !ja = k_ij + ja = 2*k_ij-1 + else + ! bottom-left corner + !ja = k_ij !FIXED BUG + ja = 2*k_ij-1 !FIXED BUG + if (bc_xlo==3) sa = -sa ! reflect U in x + endif + + if (ja > -1) then + ! add a new matrix element for interior, wall, Neumann case + nelt = nelt+1 + minfo_matrix_ia(nelt) = 2*k_ij + minfo_matrix_ja(nelt) = ja + minfo_matrix_sa(nelt) = sa + endif + + + !------------------------------------------------------- + ! Right hand sides for vector elements k_ij and k_ij+nD + ! Need to add comments on SGN RHS. + + if (debug) then + write(67,*) '++before i,j,rhs = ',i,j,rhs(2*k_ij-1),rhs(2*k_ij) + endif + + !RHS for SGN + phix = (phi(i+1,j) - phi(i-1,j))/(2.d0*dxm) + phiy = (phi(i,j+1) - phi(i,j-1))/(2.d0*dym) + + wx = (w(i+1,j) - w(i-1,j))/(2.d0*dxm) + wy = (w(i,j+1) - w(i,j-1))/(2.d0*dym) + + etax = (eta(i+1,j) - eta(i-1,j))/(2.d0*dxm) + etay = (eta(i,j+1) - eta(i,j-1))/(2.d0*dym) + + rhs(2*k_ij-1) = rhs(2*k_ij-1) + (grav/alpha * etax & + + 2.d0*h_ij*(h_ij/3.d0 * phix & + + phi(i,j) * (hx + 0.5d0*Bx(i,j))) & + + 0.5d0*h_ij*wx + w(i,j)*etax) + + rhs(2*k_ij) = rhs(2*k_ij) + (grav/alpha * etay & + + 2.d0*h_ij*(h_ij/3.d0 * phiy & + + phi(i,j) * (hy + 0.5d0*By(i,j))) & + + 0.5d0*h_ij*wy + w(i,j)*etay) + + if (debug) then + write(67,*) '++after i,j,rhs = ',i,j,rhs(2*k_ij-1),rhs(2*k_ij) + endif + + if (debug) then + write(44,998) nelt,nelt-neltSt + 998 format(" ending with nelt ",i6," used ",i5) + endif + + + if (debug) then ! dump matrix to look for singularity + write(88,*)" triplets for i,j ",i,j," level ",levelBouss + do k = neltBegin+1, nelt + !write(88,103) minfo_matrix_ia(k),minfo_matrix_ja(k),minfo_matrix_sa(k) + write(88,103) minfo_matrix_sa(k) +! 103 format(2i8,e16.7) + 103 format(e16.7) + end do + write(88,*) + !close(88) + !write(89,*)" level ",levelBouss + !do k = 1,2*numBoussCells + ! write(89,104) k,rhs(k) + 104 format(i5,e16.7) + !end do + !!close(89) + !!stop + endif + + enddo + enddo + + if (debug) then + if (mptr .eq. 1) then + write(88,*)" grid ",mptr," level ",levelBouss,"x,y,ia,ja,sa" + !do ii = 0, nx+1 + !do jj = 0, ny+1 + do ii = 1, nx + do jj = 1, ny + k = mi%mindex(ii,jj) + kd = k + nD + x = xlow + (ii-.5d0)*dx + y = ylow + (jj-.5d0)*dy + write(88,105) x,y,k,minfo_matrix_ia(k),minfo_matrix_ja(k),minfo_matrix_sa(k) + write(88,105) x,y,k+nD,minfo_matrix_ia(k+nD),minfo_matrix_ja(k+nD),minfo_matrix_sa(k+nD) + 105 format(2e15.7,3i8,e16.7) + 799 format(2i5,4e15.7) + 798 format(6e15.7,3i5) + end do + end do + endif + endif + + + + ! put back in struct. May have more grids following + minfo%matrix_nelt = nelt + +#ifdef WHERE_AM_I + write(*,*) "ending buildSparseMatrixSGNcoo" +#endif + + + return +end subroutine buildSparseMatrixSGNcoo_blocks diff --git a/src/2d/bouss/implicit_update_bouss_2Calls.f90 b/src/2d/bouss/implicit_update_bouss_2Calls.f90 index b0af86ffc..3112dd6b0 100644 --- a/src/2d/bouss/implicit_update_bouss_2Calls.f90 +++ b/src/2d/bouss/implicit_update_bouss_2Calls.f90 @@ -10,7 +10,7 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time) ! Boussinesq version. use amr_module - use topo_module, only: topo_finalized + use topo_module, only: aux_finalized use bouss_module implicit none @@ -62,10 +62,11 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time) !!write(*,*) "before bound loop" !!call timestamp() - !! reset isBouss flag array if topo still moving - !if (.not. topo_finalized) then + !! reset isBouss flag array if topo/aux still changing + + if (aux_finalized .lt. 2) then call setBoussFlag(levelBouss,naux) - !endif + endif !$OMP PARALLEL DO PRIVATE(j,locnew,locaux,mptr,nx,ny,mitot, & !$OMP mjtot,levSt), & !$OMP SHARED(levelBouss,nvar,naux,alloc,nghost,time, & @@ -127,7 +128,7 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time) ! ! convert to compressed row format ! nst = minfo%matrix_nelt ! call st_to_cc_size(nst,minfo%matrix_ia,minfo%matrix_ja,ncc) - ! call pardiso_driver(soln,rhs,levelBouss,numBoussCells,nst,ncc,topo_finalized) + ! call pardiso_driver(soln,rhs,levelBouss,numBoussCells,nst,ncc,aux_finalized) ! call system_clock(clock_finishLinSolve,clock_rate) ! call cpu_time(cpu_finishLinSolve) ! timeLinSolve = timeLinSolve + clock_finishLinSolve - clock_startLinSolve @@ -144,7 +145,7 @@ subroutine implicit_update(nvar,naux,levelBouss,numBoussCells,doUpdate,time) call system_clock(clock_startLinSolve,clock_rate) call cpu_time(cpu_startLinSolve) call petsc_driver(soln,rhs,levelBouss,numBoussCells,time, & - topo_finalized,minfo%matrix_nelt) + aux_finalized,minfo%matrix_nelt) call system_clock(clock_finishLinSolve,clock_rate) call cpu_time(cpu_finishLinSolve) !write(89,*)" level ",levelBouss," rhs soln" @@ -416,10 +417,10 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & use amr_module, only: outunit,mxnest use bouss_module, only: matrix_info_allLevs,boussMindepth,crs use bouss_module, only: matrix_levInfo, matrix_patchIndex + use bouss_module, only: alpha, origCooFormat use amr_module, only: rnode,cornxlo,cornylo,hxposs,hyposs use geoclaw_module, only: earth_radius, deg2rad, coordinate_system, grav use geoclaw_module, only: dry_tolerance - use bouss_module, only: alpha implicit none integer, intent(in) :: nvar,naux,nghost,mitot,mjtot,nD,mptr,nb,levelBouss @@ -466,10 +467,13 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & !debug = .true. debug = .false. - + !if (levelBouss .eq. 1) debug = .true. !! if not updating momentum, just save psi, translating from !! soln vector back to patch variables + if (debug) then + write(*,*)"start of solnUpdateSGN, level ",levelBouss," doUpate = ",doUpdate + endif !! if finest level grid, dont have 2 levels of storage !! no need to save in old since wont need for interpolation @@ -506,18 +510,35 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & if (mi%isBouss(i-nghost,j-nghost)) then ! check if soln == 0 as proxy if cell is more shallow than deep bouss ! and no update - if (.not. crs) then + if (.not. crs) then ! COO triplet format - valnew(2,i,j) = valnew(2,i,j) + dt * valnew(1,i,j)* & - (grav/alpha * etax(i,j) - soln(k_ij)) - valnew(3,i,j) = valnew(3,i,j) + dt * valnew(1,i,j)* & - (grav/alpha * etay(i,j) - soln(k_ij+nD)) - else ! CRS format, different mapping + if (origCooFormat) then + valnew(2,i,j) = valnew(2,i,j) + dt * valnew(1,i,j)* & + (grav/alpha * etax(i,j) - soln(k_ij)) + valnew(3,i,j) = valnew(3,i,j) + dt * valnew(1,i,j)* & + (grav/alpha * etay(i,j) - soln(k_ij+nD)) + if (debug) then + write(*,901) i,j,k_ij,soln(k_ij),soln(k_ij+nD),rhs(k_ij),rhs(k_ij+nD) + endif + else ! block format, 1-based + 901 format(3i9,4e18.8) + valnew(2,i,j) = valnew(2,i,j) + dt * valnew(1,i,j)* & + (grav/alpha * etax(i,j) - soln(2*k_ij-1)) + valnew(3,i,j) = valnew(3,i,j) + dt * valnew(1,i,j)* & + (grav/alpha * etay(i,j) - soln(2*k_ij)) + if (debug) then + write(*,901) i,j,k_ij,soln(2*k_ij-1),soln(2*k_ij),rhs(2*k_ij-1),rhs(2*k_ij) + endif + endif + else ! CRS format, different mapping, 0 based valnew(2,i,j) = valnew(2,i,j) + dt * valnew(1,i,j)* & (grav/alpha * etax(i,j) - soln(2*(k_ij-1))) valnew(3,i,j) = valnew(3,i,j) + dt * valnew(1,i,j)* & (grav/alpha * etay(i,j) - soln(2*k_ij-1)) - endif + if (debug) then + write(*,901)i,j,k_ij,soln(2*(k_ij-1)),soln(2*k_ij-1),rhs(2*(k_ij-1)),rhs(2*k_ij-1) + endif + endif ! save update in components 4,5 for Dirichlet BC next iteration ! for coarser levels, valOther is valOld @@ -525,9 +546,14 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & ! and valOther is valNew, since dont need to save old ! and new for interpolation if (.not. crs) then - valOther(4,i,j) = soln(k_ij) - valOther(5,i,j) = soln(k_ij+nD) - else + if (origCooFormat) then + valOther(4,i,j) = soln(k_ij) + valOther(5,i,j) = soln(k_ij+nD) + else ! block format + valOther(4,i,j) = soln(2*k_ij-1) + valOther(5,i,j) = soln(2*k_ij) + endif + else ! crs format !valOther(4,i,j) = soln(2*k_ij-1) !valOther(5,i,j) = soln(2*k_ij) valOther(4,i,j) = soln(2*(k_ij-1)) @@ -547,6 +573,7 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & ! ghost cells should be updated. CORNERS? ! remember mindex only from 1:nx, 1:ny so subtract nghost as needed + ! below code left in to remember why it doesnt work go to 444 ! left side do j = nghost+1, mjtot-nghost @@ -616,21 +643,35 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & k_ij = mi%mindex(i-nghost,j-nghost) ! hu element index for (i,j) if (mi%isBouss(i-nghost,j-nghost)) then if (.not. crs) then - valnew(4,i,j) = soln(k_ij) - valnew(5,i,j) = soln(k_ij+nD) - else - valnew(4,i,j) = soln(2*(k_ij-1)) - valnew(5,i,j) = soln(2*k_ij-1) + if (origCooFormat) then + valnew(4,i,j) = soln(k_ij) + valnew(5,i,j) = soln(k_ij+nD) + if (debug) then + write(*,901)i,j,k_ij,soln(k_ij),soln(k_ij+nD),rhs(k_ij),rhs(k_ij+nD) + endif + else ! block format + valnew(4,i,j) = soln(2*k_ij-1) + valnew(5,i,j) = soln(2*k_ij) + if (debug) then + write(*,901)i,j,k_ij,soln(2*k_ij-1),soln(2*k_ij),rhs(2*k_ij-1),rhs(2*k_ij) + endif + endif + else ! using crs + valnew(4,i,j) = soln(2*(k_ij-1)) + valnew(5,i,j) = soln(2*k_ij-1) + if (debug) then + write(*,901)i,j,k_ij,soln(2*(k_ij-1)),soln(2*k_ij-1),rhs(2*(k_ij-1)),rhs(2*k_ij-1) + endif endif - else + else ! not bouss cell valnew(4,i,j) = 0.d0 valnew(5,i,j) = 0.d0 endif - end do - end do + end do + end do endif - if (debug .or. 0 .eq. 1) then + if (debug) then xlowWithGhost = rnode(cornxlo,mptr) - nghost*dx ylowWithGhost = rnode(cornylo,mptr) - nghost*dy write(34,*)" end of solnUpdateSGN for grid ",mptr," with doUpdate ",doUpdate @@ -640,29 +681,6 @@ subroutine solnUpdateSGN(valnew,valOther,nb,mitot,mjtot,nghost,nvar,dt, & mitot,mjtot,ng,iunit,doUpdate,levelBouss) endif - if (.true. .and. debug) then - imax = -1 - jmax = -1 - speedmax = -1.d0 - do j = nghost+1, mjtot-nghost - do i = nghost+1, mitot-nghost - if (valnew(1,i,j) .gt. dry_tolerance) then - speed(i,j) = sqrt((valnew(2,i,j)**2+valnew(3,i,j)**2))/valnew(1,i,j) - if (speed(i,j) .gt. speedmax) then - speedmax = speed(i,j) - imax = i - jmax = j - endif - else - speed(i,j) = 0.d0 - endif - end do - end do - write(*,333) mptr,mitot,mjtot,speedmax,imax,jmax, & - (valnew(ivar,imax,jmax),ivar=1,nvar) - 333 format("grid ",i4," dims ",i4,i4," speed",f9.2," max index at ",2i4,5e12.4) - endif - end subroutine solnUpdateSGN subroutine symcheck(val,mitot,mjtot) diff --git a/src/2d/bouss/petsc_driver.f90 b/src/2d/bouss/petsc_driver.f90 index 12be87372..b1ac666cd 100644 --- a/src/2d/bouss/petsc_driver.f90 +++ b/src/2d/bouss/petsc_driver.f90 @@ -1,4 +1,5 @@ -subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalized,nelt) +subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time, & + aux_finalized,nelt) ! use petsc to solve sparse linear system ! called this way so can allocate storage of correct size @@ -13,8 +14,7 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize integer, intent(in) :: levelBouss, numBoussCells real(kind=8), intent(inout) :: rhs_geo(0:2*numBoussCells) real(kind=8), intent(in) :: time - logical, intent(in) :: topo_finalized - integer, intent(in) :: nelt + integer, intent(in) :: nelt, aux_finalized !! pass in numBoussCells for this level so can dimension this array real(kind=8), intent(out) :: soln(0:2*numBoussCells) @@ -56,7 +56,7 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize !write(17,777)(j+1,minfo%matrix_ia(j)+1,minfo%matrix_ja(j)+1,minfo%matrix_sa(j),j=0,nelt-1) 777 format(3i8,e15.6) - if (newGrids(levelBouss) .or. .not. topo_finalized) then + if (newGrids(levelBouss) .or. aux_finalized .lt. 2) then ! destroy old solver objects to recreate for new grid setup. No need ! to check if init or restart because petsc checks if null call KSPDestroy(ksp(levelBouss),ierr) @@ -168,7 +168,7 @@ subroutine petsc_driver(soln,rhs_geo,levelBouss,numBoussCells,time,topo_finalize numTimes(levelBouss) = numTimes(levelBouss) + 1 !DEBUG - !call vecDuplicate(solution,y,ierr) !makes new vec y same size as first arg + !call vecDuplicate(solution,y,ierr) ! makes new vec y with same size as soln !call MatMult(Jr(levelBouss),solution,y,ierr) !call VecView(rhs,PETSC_VIEWER_STDOUT_SELF,ierr) !call VecView(y,PETSC_VIEWER_STDOUT_SELF,ierr) diff --git a/src/2d/bouss/prepBuildSparseMatrixSGNcoo.f90 b/src/2d/bouss/prepBuildSparseMatrixSGNcoo.f90 index 35f7154c5..dd120f058 100644 --- a/src/2d/bouss/prepBuildSparseMatrixSGNcoo.f90 +++ b/src/2d/bouss/prepBuildSparseMatrixSGNcoo.f90 @@ -52,7 +52,7 @@ subroutine prepBuildSparseMatrixSGNcoo(soln,rhs,nvar,naux,levelBouss,numBoussCel call system_clock(clock_start,clock_rate) call cpu_time(cpu_start) - !debug = .false. + debug = .false. !debug = .true. minfo => matrix_info_allLevs(levelBouss) @@ -101,24 +101,37 @@ subroutine prepBuildSparseMatrixSGNcoo(soln,rhs,nvar,naux,levelBouss,numBoussCel if (locold .eq. 0) locold = loc ! finest level has no old time soln locaux = node(storeaux,mptr) ! needed for iaddaux - call buildSparseMatrixSGNcoo(alloc(loc),alloc(locold), & + if (origCooFormat) then + call buildSparseMatrixSGNcoo(alloc(loc),alloc(locold), & alloc(locaux), soln,rhs, & minfo%matrix_ia,minfo%matrix_ja,minfo%matrix_sa, & numBoussCells,levelBouss,mptr,nx,ny,nvar,naux, & matrix_rhs1_max,matrix_rhs2_max) + else + call buildSparseMatrixSGNcoo_blocks(alloc(loc),alloc(locold), & + alloc(locaux), soln,rhs, & + minfo%matrix_ia,minfo%matrix_ja,minfo%matrix_sa, & + numBoussCells,levelBouss,mptr,nx,ny,nvar,naux, & + matrix_rhs1_max,matrix_rhs2_max) + endif enddo nng_loop - !if (minfo%matrix_nelt .gt. 0) then - !! copy into allocatable to continue. stack vars needed for debugging - ! minfo%matrix_ia = minfo_matrix_ia - ! minfo%matrix_ja = minfo_matrix_ja - ! minfo%matrix_sa = minfo_matrix_sa - ! !write(35,909) (j,minfo_matrix_ia(j),minfo_matrix_ja(j),minfo_matrix_sa(j),j=1,minfo%matrix_nelt) - 909 format(3i6,e25.15) - !endif + if (minfo%matrix_nelt .gt. 0 .and. debug) then + ! output for testing in matlab + !minfo%matrix_ia = minfo_matrix_ia + !minfo%matrix_ja = minfo_matrix_ja + !minfo%matrix_sa = minfo_matrix_sa + write(25,*)" coo matrix level ",levelBouss," num rows ",minfo%matrix_nelt + write(25,909) (j,minfo%matrix_ia(j),minfo%matrix_ja(j),minfo%matrix_sa(j),j=1,minfo%matrix_nelt) + 909 format(3i8,e25.15) + write(26,*)" rhs size",2*numBoussCells," u and v equations" + write(26,908) (j,rhs(j),j=1,numBoussCells) + write(26,908) (j+nD,rhs(j+nD),j=1,numBoussCells) + 908 format(i8,e25.15) + endif !============ Debug: diff --git a/src/python/geoclaw/fgout_tools.py b/src/python/geoclaw/fgout_tools.py index 862eeeacc..2f6849b53 100644 --- a/src/python/geoclaw/fgout_tools.py +++ b/src/python/geoclaw/fgout_tools.py @@ -100,7 +100,7 @@ def qmap(self): @property def h(self): - """depth""" + """Depth""" if self._h is None: q_out_vars = self.fgout_grid.q_out_vars try: @@ -118,7 +118,7 @@ def h(self): @property def hu(self): - """momentum h*u""" + """Momentum h*u""" if self._hu is None: q_out_vars = self.fgout_grid.q_out_vars try: @@ -131,7 +131,7 @@ def hu(self): @property def u(self): - """speed u, computed as hu/h or set to 0 if h 1: + ax.add_patch(plt.Circle((x[-1], y[-1]), _radius[-1], + color=fill_color)) + + # Draw path around inner points + if t.shape[0] > 2: + for i in range(t.shape[0] - 1): + p = np.array([(x[i], y[i]), (x[i + 1], y[i + 1])]) + v = p[1] - p[0] + if abs(v[1]) > 1e-16: + n = np.array([1, -v[0] / v[1]], dtype=float) + elif abs(v[0]) > 1e-16: + n = np.array([-v[1] / v[0], 1], dtype=float) + else: + continue + # raise Exception("Zero-vector given") + n /= np.linalg.norm(n) + n *= _radius[i] + + ax.fill((p[0, 0] + n[0], p[0, 0] - n[0], + p[1, 0] - n[0], + p[1, 0] + n[0]), + (p[0, 1] + n[1], p[0, 1] - n[1], + p[1, 1] - n[1], + p[1, 1] + n[1]), + facecolor=fill_color, alpha=fill_alpha) + ax.add_patch(plt.Circle((p[1][0], p[1, 1]), _radius[i], + color=fill_color, alpha=fill_alpha)) + + def write_hwrf(self, path, verbose=False): + r"""Write out an TCVITALS formatted storm file + + :Input: + - *path* (string) Path to the file to be written. + - *verbose* (bool) Print out additional information when writing. + """ + raise NotImplementedError("HWRF formatted info files cannot be ", + "written out yet.") + + def write_owi(self, path, verbose=False): + r"""Write out an OWI information formatted storm file + + :Input: + - *path* (string) Path to the file to be written. + - *verbose* (bool) Print out additional information when writing. + """ + # OWI file formats + _data_file_format_mapping = {'ascii': 1, 'nws12': 1, + 'netcdf': 2, 'nws13': 2} + + try: + with path.open('w') as data_file: + # Write header + data_file.write("# OWI Data Decription - NWS12 and NWS13\n") + if isinstance(self.data_file_format, int): + file_format = self.data_file_format + elif isinstance(self.data_file_format, str): + if (self.data_file_format.lower() in + _data_file_format_mapping.keys()): + file_format = _data_file_format_mapping[ + self.data_file_format.lower()] + else: + raise TypeError("Unknown storm data file format type" + + f" '{self.data_file_format}' provided.") + else: + raise TypeError("Unknown storm data file format type" + + f" '{self.data_file_format}' provided.") + + # Write out data + # Time offset + self.time_offset = np.datetime64(self.time_offset) + if isinstance(self.time_offset, np.datetime64): + t = np.datetime_as_string(self.time_offset, unit="s") + data_file.write(f"{t.ljust(20)} # Time Offset\n") + else: + raise ValueError("Time offset must be a datetime64 object.") + + # File format + data_file.write(f"{str(file_format).ljust(20)} # File format\n") + data_file.write("\n") + + # File format specific values + data_file.write("# Data Information\n") + if file_format == 1: + # :TODO: Modify number of regions to be indpendent of file_paths + if len(self.file_paths)%2 != 0: + raise ValueError("The number of files should be even, " + + "one for pressure and wind, for each " + + "resolution provided.") + num_regions = int(len(self.file_paths) / 2) + data_file.write(f"{str(num_regions).ljust(20)}" + + " # Number of regions\n") + for n in range(int(len(self.file_paths) / 2)): + data_file.write(f"{str(self.file_paths[n * 2]).ljust(20)}\n") + data_file.write(f"{str(self.file_paths[n * 2 + 1]).ljust(20)}\n") + elif file_format == 2: + if len(self.file_paths) != 1: + raise ValueError("Expected 1 file for format " + + f"{file_format} rather than " + + f"{len(self.file_paths)}") + data_file.write(f"{self.file_paths[0]}\n") + + except Exception as e: + # If an exception occurs clean up a partially generated file + Path.unlink(path, missing_ok=True) + raise e + + # ================ + # Track Plotting + # ================ + def plot(self, ax, *args, t_range=None, categorization=None, + cat_colors={}, plot_swath=False, radius=None, + coordinate_system=2, fill_alpha=0.25, fill_color='red', + **kwargs): + """Plot this storm's track in the given axes object + + :Input: + - *ax* (matplotlib.pyplot.axes) Axes to plot into. + - *t_range* (list) Time range to plot the track for. If None then use + entire range. Default is None. + - *categorization* (str) Type of categorization to be used. This is + used to map to the keys in the cat_colors dictionary. Default is + None and will cause no categorization to occur. + - *cat_colors* (dict) Color mapping between numeric categorization and + colors to be plotted for the track. + - *plot_swath* (bool) Plot a swath around the track using one of the + methods determined by what radius information is provided. Default + is False. + - *radius* (None or float or numpy.ndarray) + - *coordinate_system* (int) + - *fill_alpha* (float) + - *fill_color* (color) + - *kwargs* All additional keyword arguments are passed to the plotting + command for the track. + """ import matplotlib.pyplot as plt # Extract information for plotting the track/swath @@ -1830,7 +2011,6 @@ def category(self, categorization="NHC", cat_names=False): *string*. This is only returned if *car_names = True*. """ - # TODO: Need to standardize on 1-minute (almost never available) or # 10-minute (widely available) - see # https://en.wikipedia.org/wiki/Tropical_cyclone#Major_basins_and_related_warning_centers @@ -2008,7 +2188,6 @@ def fill_rad_w_other_source(t, storm_targ, storm_fill, var, interp_kwargs={}): ... path = 'out_path.storm', ... max_wind_radius_fill = fill_mwr) """ - fill_da = xr.DataArray(getattr(storm_fill, var), coords={'t': getattr(storm_fill, 't')}, dims=('t',)) @@ -2204,6 +2383,7 @@ def write_data_derived(self, filename): :param filename: The name of the output NetCDF file (without the extension). :return: None """ + import xarray as xr windx = numpy.array(self.u) windy = numpy.array(self.v)