Transport: added an option to specify spatially varying epsilon2 for …

…weno.

sample_inputs/param.nml

@@ -515,13 +515,23 @@
   courant_weno=0.5 !Courant number for weno transport
   nquad = 2  !number of quad points on each side, nquad= 1 or 2
   ntd_weno = 1 !order of temporal discretization: (1) Euler (default); (3): 3rd-order Runge-Kutta (only for benchmarking)
-  epsilon1 = 1.e-15   !coefficient for 2nd order weno smoother (larger values are more prone to numerical dispersion)
-  epsilon2 = 1.e-10  !1st coefficient for 3rd order weno smoother (larger values are more prone to numerical dispersion
-                     !, 1.e-10 should be fairly safe, recommended values: 1.e-8 ~ 1.e-6 (from the real applications so far)
+
+  epsilon1 = 1.e-15   !coefficient of 2nd-order weno smoother (larger values are more prone to numerical dispersion)
+  i_epsilon2 = 1      !switch for the specification method of 3rd-order weno smoother coefficient
+                      !(1) spatially uniform, set by the parameter epsilon2 below;
+                      !(2) spatially varying, set by the input file epsilon2.gr3,
+                      !    Important:
+                      !      Provide log10(epsilon2) in epsilon2.gr3, e.g.,
+                      !      if epsilon2=10^(-6) or 1e-6, then write the z value should be -6 in epsilon2.gr3.
+                      !      This is mainly for the convenience of visualization.
+  epsilon2 = 1.e-10   !used when i_epsilon2 = 1
+                      !spatially uniform coefficient of 3rd-order weno smoother
+                      !(larger values are more prone to numerical dispersion,
+                      ! 1.e-10 should be fairly safe, recommended values: 1.e-8 ~ 1.e-6)
   i_prtnftl_weno = 0 !option for writing nonfatal errors on invalid temp. or salinity for density: (0) off; (1) on.
 
   !Inactive at the moment:
-  epsilon3 = 1.e-25  !2nd coefficient for 3rd order weno smoother (inactive at the moment)
+  epsilon3 = 1.e-25  !Second coefficient of 3rd-order weno smoother (inactive at the moment)
   !Elad filter has not been implemented yet; preliminary tests showed it might not be necessary
   ielad_weno = 0      !ielad, if ielad=1, use ELAD method to suppress dispersion (inactive at the moment)
   small_elad = 1.e-4  !small (inactive at the moment)

src/Core/schism_glbl.F90

@@ -496,17 +496,19 @@ module schism_glbl
   integer,save,allocatable :: isbe(:) !(ne): bnd seg flags, isbe(ie)=1 if any node of element ie lies on bnd; isbe(ie)=0 otherwise
   logical,save,allocatable :: is_inter(:)  !identifier of interface sides (between two ranks), for debugging only
   integer,save,allocatable :: iside_table(:) !a record of all interface sides within the current rank
-  integer, save :: ip_weno !order of the polynomials used for weno stencils, see param.in.sample
+  integer,save :: ip_weno !order of the polynomials used for weno stencils, see param.in.sample
   real(rkind),save :: courant_weno !Courant number for weno transport
   real(rkind),save :: epsilon1 !coefficient for 2nd order weno smoother
-  real(rkind),save :: epsilon2 !1st coefficient for 3rd order weno smoother
+  integer,save :: i_epsilon2 !switch for 1st coefficient of 3rd order weno smoother
+  real(rkind),save :: epsilon2 !1st coefficient of 3rd order weno smoother, spatially uniform value from param.nml
+  real(rkind),save,allocatable :: epsilon2_elem(:) !1st coefficient of 3rd order weno smoother, elemental values
   real(rkind),save :: epsilon3 !2nd coefficient for 3rd order weno smoother
-  integer, save :: nquad !number of quad points used for 3rd order weno
+  integer,save :: nquad !number of quad points used for 3rd order weno
   !levels of time discretization, mainly for testing purposes
-  integer, save :: ntd_weno !(1) one-level, reduces to Euler; (2) not implemented yet; (3) 3rd-order Runge-Kutta temporal discretization (Shu and Osher, 1988)
+  integer,save :: ntd_weno !(1) one-level, reduces to Euler; (2) not implemented yet; (3) 3rd-order Runge-Kutta temporal discretization (Shu and Osher, 1988)
   !Elad filter
-  integer, save :: ielad_weno !switch for elad filter, not used at the moment
-  integer, save :: i_prtnftl_weno !switch for printing invalid T/S to nonfatal_*
+  integer,save :: ielad_weno !switch for elad filter, not used at the moment
+  integer,save :: i_prtnftl_weno !switch for printing invalid T/S to nonfatal_*
   real(rkind),save :: small_elad !criteria for ELAD, not used at the moment
 
   real(rkind),save,allocatable :: xqp(:,:),yqp(:,:)  !quadrature point coordinates

src/Hydro/schism_init.F90

@@ -196,7 +196,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
      &vdmax_pp2,vdmin_pp1,vdmin_pp2,tdmin_pp1,tdmin_pp2,mid,stab,xlsc0, &
      &ibcc_mean,flag_ic,start_year,start_month,start_day,start_hour,utc_start, &
      &itr_met,h_tvd,eps1_tvd_imp,eps2_tvd_imp,ip_weno, &
-     &courant_weno,ntd_weno,nquad,epsilon1,epsilon2,epsilon3,ielad_weno,small_elad, &
+     &courant_weno,ntd_weno,nquad,epsilon1,i_epsilon2,epsilon2,epsilon3,ielad_weno,small_elad, &
      &i_prtnftl_weno,inu_tr,step_nu_tr,vnh1,vnh2,vnf1,vnf2, &
      &moitn0,mxitn0,rtol0,iflux,inter_mom,h_bcc1,inu_elev,inu_uv, &
      &ihhat,kr_co,rmaxvel,velmin_btrack,btrack_nudge,ibtrack_test,irouse_test, &
@@ -467,7 +467,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
       mid='KL'; stab='KC'; xlsc0=0.1_rkind;  
       ibcc_mean=0; flag_ic(:)=1; start_year=2000; start_month=1; start_day=1; start_hour=0._rkind; utc_start=8._rkind;  
       itr_met=3; h_tvd=5._rkind; eps1_tvd_imp=1.d-4; eps2_tvd_imp=1.d-14; ip_weno=2;  
-      courant_weno=0.5_rkind; ntd_weno=1; nquad=2; epsilon1=1.d-15; epsilon2=1.d-10; epsilon3=1.d-25; 
+      courant_weno=0.5_rkind; ntd_weno=1; nquad=2; epsilon1=1.d-15; i_epsilon2=1; epsilon2=1.d-10; epsilon3=1.d-25; 
       ielad_weno=0; small_elad=1.d-4; i_prtnftl_weno=0;
       inu_tr(:)=0; step_nu_tr=86400._rkind; vnh1=400._rkind; vnh2=500._rkind; vnf1=0._rkind; vnf2=0._rkind;
       moitn0=50; mxitn0=1500; rtol0=1.d-12; iflux=0; inter_mom=0; 
@@ -870,12 +870,18 @@ subroutine schism_init(iorder,indir,iths,ntime)
           call parallel_abort(errmsg)
         endif
 
-!        call get_param('param.in','epsilon2',2,itmp,epsilon2,stringvalue)
-        if(epsilon2<=0._rkind) then
-          write(errmsg,*)'Illegal epsilon2:',epsilon2
+        if(i_epsilon2.ne.1 .and. i_epsilon2.ne.2) then
+          write(errmsg,*)'Illegal i_epsilon2:',i_epsilon2
           call parallel_abort(errmsg)
         endif
 
+        if (i_epsilon2.eq.1) then
+          if(epsilon2<=0._rkind) then
+            write(errmsg,*)'Illegal epsilon2:',epsilon2
+            call parallel_abort(errmsg)
+          endif
+        endif
+
 !        call get_param('param.in','epsilon3',2,itmp,epsilon3,stringvalue)
         if(epsilon3<=0._rkind) then
           write(errmsg,*)'Illegal epsilon3:',epsilon3
@@ -1412,7 +1418,7 @@ subroutine schism_init(iorder,indir,iths,ntime)
          &  sparsem(0:mnei_p,np), & !sparsem for non-ghosts only
          &  tr_nudge(natrm,npa), & 
          &  fun_lat(0:2,npa),dav(2,npa),elevmax(npa),dav_max(2,npa),dav_maxmag(npa), &
-         &  diffmax(npa),diffmin(npa),dfq1(nvrt,npa),dfq2(nvrt,npa), & 
+         &  diffmax(npa),diffmin(npa),dfq1(nvrt,npa),dfq2(nvrt,npa),epsilon2_elem(ne), & 
          &  iwater_type(npa),rho_mean(nvrt,nea),erho(nvrt,nea),& 
          & surf_t1(npa),surf_t2(npa),surf_t(npa),etaic(npa),sav_alpha(npa), &
          & sav_h(npa),sav_nv(npa),sav_di(npa),sav_cd(npa), &
@@ -3586,6 +3592,40 @@ subroutine schism_init(iorder,indir,iths,ntime)
         enddo !i
       endif
 
+!     epsilon2 (1st coefficient) of 3rd order weno smoother
+      if (i_epsilon2.eq.1) then
+          epsilon2_elem(:)=epsilon2
+      elseif (i_epsilon2.eq.2) then
+        if(myrank==0) then
+          open(32,file=in_dir(1:len_in_dir)//'epsilon2.gr3',status='old')
+          read(32,*)
+          read(32,*) itmp1,itmp2
+          if(itmp1/=ne_global.or.itmp2/=np_global) & 
+     &call parallel_abort('Check epsilon2.gr3')
+          do i=1,np_global
+            read(32,*)j,xtmp,ytmp,tmp
+            buf3(i)=tmp
+          enddo !i
+          close(32)
+        endif !myrank
+        call mpi_bcast(buf3,ns_global,rtype,0,comm,istat)
+
+        do i=1,np_global
+          if(ipgl(i)%rank==myrank) then
+            swild(ipgl(i)%id)=buf3(i)
+          endif
+        enddo !i
+
+        do i=1,ne
+          tmp = minval(swild(elnode(1:i34(i),i)))  !smaller values are safer
+          if (tmp>-2.0) then
+            write(errmsg,*)'Invalid epsilon2 (must <=0.01)', 10.0**tmp, ' at Element ', ielg(i)
+            call parallel_abort(errmsg)
+          endif
+          epsilon2_elem(i)=10.0**tmp
+        enddo !i
+      endif
+
 !...  Land b.c. option (inactive)
 !      read(15,*) !islip !0: free slip; otherwise no slip
       islip=0

src/Hydro/transport_TVD_imp.F90

@@ -581,9 +581,9 @@ subroutine do_transport_tvd_imp(it,ntr,difnum_max_l) !,nvrt1,npa1,dfh1)
                   endif
 #endif
 
-                  wm(i)=1.0d0/((epsilon2+b4)*(epsilon2+b4))
+                  wm(i)=1.0d0/((epsilon2_elem(ie)+b4)*(epsilon2_elem(ie)+b4))
                   !----test-------------- 
-                  !wm1(i)=1.0/((epsilon2+b1+b2)*(epsilon2+b1+b2)); ! sum1=sum1+wm1(i)
+                  !wm1(i)=1.0/((epsilon2_elem(ie)+b1+b2)*(epsilon2_elem(ie)+b1+b2)); ! sum1=sum1+wm1(i)
                   !wm2(i)=1.0/((epsilon3+b3)*(epsilon3+b3)); 
                   !wm(i)=wm1(i)*wm2(i)
                   !----------------------