From 4f44788d43cd14e7fba0a8bb4ad245faee6f0bdf Mon Sep 17 00:00:00 2001 From: Joseph Zhang Date: Thu, 4 Jan 2024 13:16:16 -0500 Subject: [PATCH] Add a script to convert SMS scatter to ACE/region format. Revised gen_tvd_WENO.f90 to allow multiple regions. --- src/Utility/Pre-Processing/gen_tvd_WENO.f90 | 150 +++++++++++++++----- src/Utility/SMS/scatter_2_region.f90 | 53 +++++++ 2 files changed, 169 insertions(+), 34 deletions(-) create mode 100644 src/Utility/SMS/scatter_2_region.f90 diff --git a/src/Utility/Pre-Processing/gen_tvd_WENO.f90 b/src/Utility/Pre-Processing/gen_tvd_WENO.f90 index 91cdfac27..ea69373e7 100644 --- a/src/Utility/Pre-Processing/gen_tvd_WENO.f90 +++ b/src/Utility/Pre-Processing/gen_tvd_WENO.f90 @@ -1,50 +1,103 @@ ! Generate tvd.prop for hybrid WENO/ELM and cross-scale applications, based on depth and remove ! 'isolated' WENO elements in shallow water, etc. The goal is to -! minimize the direct 'contact' between WENO and ELM cells, which -! can cause dispersion. +! minimize the direct 'contact' between WENO and ELM cells, which can cause dispersion. ! Inputs: (1) screen; ! (2) hgrid.gr3 (in any projection or lon/lat; b.c. part not needed) -! (3) nearshore.gr3: '0' for eddying regime; non-0 for nearshore. Nearshore region -! should be slightly offshore of isobath 'hmin_of' -! Output: tvd.prop.0 (may be further edited, e.g. deep channels/fjords) +! (3) nearshore.rgn: nearshore region (e.g. using 20m isoth) inside +! which upwind is used except in estuary_*.rgn; +! (4) optional estuary_[1,2..].rgn: specific estuary regions where TVD/WENO is desired; +! make sure these verlap with nearshore.rgn (so no gap in TVD/WENO zone) +! Output: tvd.prop.0 (may be further edited, e.g. connectivity etc) -! ifort -O2 -o gen_tvd_WENO ../UtilLib/schism_geometry.f90 gen_tvd_WENO.f90 +! ifort -O2 -o gen_tvd_WENO ../UtilLib/schism_geometry.f90 ../UtilLib/pt_in_poly_test.f90 gen_tvd_WENO.f90 use schism_geometry_mod + use pt_in_poly_test implicit real*8(a-h,o-z) + real*8, parameter :: small1=1.d-4 + real*8, parameter :: ray_angle=3.13192d0 + + character(len=100) :: it_char integer :: nwild(3) integer, allocatable :: i34(:),elnode(:,:),nne(:),indel(:,:),nnp(:), & - &indnd(:,:),isbnd(:),itvd(:),itvd0(:),iest(:),iest_e(:) - integer, allocatable :: ic3(:,:),elside(:,:),isdel(:,:),isidenode(:,:) + &indnd(:,:),isbnd(:),itvd(:),itvd0(:),iest(:),iest_e(:),inear(:),inear_e(:) + integer, allocatable :: ic3(:,:),elside(:,:),isdel(:,:),isidenode(:,:),nvertices(:) real*8, allocatable :: xnd(:),ynd(:),dp(:),area(:),xpoly(:),ypoly(:),xcj(:),ycj(:) - print*, 'Input cut-off depth (meters) for offshore (e.g., 30m));' - print*, 'Upwind will be used shallower than this depth in offshore:' - read*, hmin_of - - print*, 'Input cut-off depth (meters) for nearshore.gr3 (usually null() + end type + type(poly_vertex),allocatable :: poly(:) + + +! print*, 'Input cut-off depth (meters) for offshore (e.g., 30m));' +! print*, 'Upwind will be used shallower than this depth in offshore:' +! read*, hmin_of + + print*, 'Input # of estuary regions (can be 0):' + read*, nestuaries + if(nestuaries<0) stop 'nestuaries<0' + + print*, 'Input cut-off depth (meters) for estuaries (e.g. h_tvd):' + read*, hmin_estu + + !Read in nearshore.rgn and estuary regions + allocate(nvertices(nestuaries+1),poly(nestuaries+1)) !first is the nearshore.rgn + open(10,file='nearshore.rgn',status='old') + read(10,*); read(10,*)npoly + if(npoly/=1) stop 'can only have 1 poly in .rgn' + read(10,*)nvertices(1) + allocate(poly(1)%xy(nvertices(1)+1,2)) !extra pt to close the region + poly(1)%xy(nvertices(1)+1,:)=0. !init + do i=1,nvertices(1) + read(10,*)poly(1)%xy(i,1:2) !xpoly(i,1),ypoly(i,1) + enddo !i + close(10) + + !Check if the region is closed + rl=(poly(1)%xy(1,1)-poly(1)%xy(nvertices(1),1))**2+(poly(1)%xy(1,2)-poly(1)%xy(nvertices(1),2))**2 + if(rl>1.d-14) then + nvertices(1)=nvertices(1)+1 + poly(1)%xy(nvertices(1),:)=poly(1)%xy(1,:) + endif + + do irgn=1,nestuaries + write(it_char,*) irgn + it_char=adjustl(it_char); itmp=len_trim(it_char) + open(10,file='estuary_'//it_char(1:itmp)//'.rgn',status='old') + read(10,*); read(10,*)npoly + if(npoly/=1) then + print*, 'Can only have 1 poly in .rgn:',irgn + stop + endif + read(10,*)nvertices(irgn+1) + allocate(poly(irgn+1)%xy(nvertices(irgn+1)+1,2)) + poly(irgn+1)%xy(nvertices(irgn+1)+1,:)=0. + do i=1,nvertices(irgn+1) + read(10,*)poly(irgn+1)%xy(i,1:2) + enddo !i + close(10) + + !Check if the region is closed + rl=(poly(irgn+1)%xy(1,1)-poly(irgn+1)%xy(nvertices(irgn+1),1))**2+ & + &+(poly(irgn+1)%xy(1,2)-poly(irgn+1)%xy(nvertices(irgn+1),2))**2 + if(rl>1.d-14) then + nvertices(irgn+1)=nvertices(irgn+1)+1 + poly(irgn+1)%xy(nvertices(irgn+1),:)=poly(irgn+1)%xy(1,:) + endif + enddo !irgn open(14,file='hgrid.gr3',status='old') - open(13,file='nearshore.gr3',status='old') read(14,*); read(14,*) ne,np - read(13,*); read(13,*) allocate(xnd(np),ynd(np),dp(np),area(ne),i34(ne),elnode(4,ne),nne(np), & - &isbnd(np),itvd(ne),itvd0(ne),iest(np),iest_e(ne)) + &isbnd(np),itvd(ne),itvd0(ne),iest(np),iest_e(ne),inear(np),inear_e(ne)) do i=1,np read(14,*) j,xnd(i),ynd(i),dp(i) - read(13,*) j,tmp,tmp,tmp2 - iest(i)=nint(tmp2) enddo !i do i=1,ne read(14,*) j,i34(i),elnode(1:i34(i),i) - - iest_e(i)=maxval(iest(elnode(1:i34(i),i))) - n1=elnode(1,i) n2=elnode(2,i) n3=elnode(3,i) @@ -96,21 +149,43 @@ ! inputs (xnd,ynd) (x,y coordinates of each node) call schism_geometry_double(np,ne,ns,xnd,ynd,i34,elnode,ic3,elside,isdel,isidenode,xcj,ycj) +! Mark nearshore & estuary nodes + iest=0 !init as outside any estuary of interest + do i=1,np + call pt_in_poly_ray_method_double(nvertices(1),small1,ray_angle, & + &poly(1)%xy(1:nvertices(1),1),poly(1)%xy(1:nvertices(1),2),xnd(i),ynd(i),in_out,inters,npoly) + if(in_out==1) then + inear(i)=1 !nearshore + else + inear(i)=0 !offshore + endif + + !Estuaries + do irgn=1,nestuaries + call pt_in_poly_ray_method_double(nvertices(irgn+1),small1,ray_angle, & + &poly(irgn+1)%xy(1:nvertices(irgn+1),1),poly(irgn+1)%xy(1:nvertices(irgn+1),2), & + &xnd(i),ynd(i),in_out,inters,npoly) + + if(in_out==1) iest(i)=1 !inside any estuaries + enddo !irgn + enddo !i=1,np + + do i=1,ne + inear_e(i)=maxval(inear(elnode(1:i34(i),i))) + iest_e(i)=minval(iest(elnode(1:i34(i),i))) + enddo !i + ! Set TVD flag itvd(:)=0 !init do i=1,ne hmin2=minval(dp(elnode(1:i34(i),i))) - if(iest_e(i)==0) then !offshore - if(hmin2>=hmin_of) itvd(i)=1 - else !nearshore - if(hmin2>=hmin_near) itvd(i)=1 - endif !iest_e + if((inear_e(i)==0.or.iest_e(i)/=0).and.hmin2>=hmin_estu) itvd(i)=1 enddo !i ! Augment upwind zone by 1 extra layer for offshore itvd0=itvd do i=1,ne - if(iest_e(i)==0.and.itvd0(i)==0) then + if(inear_e(i)==0.and.itvd0(i)==0) then do j=1,i34(i) nd=elnode(j,i) do m=1,nne(nd) @@ -121,14 +196,14 @@ endif !ifl enddo !i -! Remove 'isolated' WENO elem nearshore +! Remove 'isolated' WENO elem loop1: do itouched=0 !# of elem changed in this iteration itvd0=itvd loop2: do i=1,ne - if(iest_e(i)==0.or.itvd0(i)==0) cycle + if(itvd0(i)==0) cycle - !WENO elem nearshore + !WENO elem icount=0 do j=1,i34(i) ie=ic3(j,i) @@ -152,9 +227,16 @@ do i=1,ne write(12,*)i,itvd(i) enddo !i - write(12,*)'hmin=',hmin_of,hmin_near + write(12,*)'hmin=',hmin_estu close(12) + ! deallocate pointer arrays + do m=1,size(poly) + if (associated(poly(m)%xy)) deallocate(poly(m)%xy) + nullify(poly(m)%xy) + enddo + deallocate(poly) + stop end diff --git a/src/Utility/SMS/scatter_2_region.f90 b/src/Utility/SMS/scatter_2_region.f90 new file mode 100644 index 000000000..443bb3227 --- /dev/null +++ b/src/Utility/SMS/scatter_2_region.f90 @@ -0,0 +1,53 @@ +! Convert SMS scatter formar .xy to .rgn (gredit), assuming it's single +! connected and in proper order +! To generate .xy in SMS: create a map and polygon (mesh generator), and edit, clean/build. Then convert +! map to Scatter (Arc end points and vertices), and highlight Scatter set and save as +! 'Scatter pt .xy'. + +! Inputs: scatter.xy (SMS) +! Outputs: scatter.rgn + +! ifort -O2 -mcmodel=medium -CB -Bstatic -o scatter_2_region scatter_2_region.f90 + + implicit real*8(a-h,o-z) + character(len=100) :: str1,str2 + real*8, allocatable :: xy(:,:) + + open(10,file='scatter.xy',status='old') + open(12,file='scatter.rgn',status='replace') + line=0 + line_start=-1 + do + read(10,*,err=99)str1 + line=line+1 + str2=adjustl(str1) +! len1=len_trim(str2) + if(str2(1:3).eq."IXY") then + line_start=line + exit + endif + enddo + +99 if(line_start==-1) stop 'did not find starting line' + print*, 'Starting line found:',line_start + + write(12,*)'Region written by ACE/gredit' + write(12,*)1 + + rewind(10) + do i=1,line_start-1; read(10,*); enddo; + read(10,*)str1,npts + write(12,*)npts,1 + allocate(xy(2,npts)) + print*, '# of points in scatter=',npts + do i=1,npts + read(10,*)j,xy(1:2,i) + write(12,*)xy(1:2,i) + enddo !i + + !Repeat 1st pt to close region: no need to do this as .rgn can be either way +! write(12,*)xy(1:2,1) + close(12) + + stop + end