TALK=T;RUN(1,1)
DISPLAY
In-Form is used to supply, at the boundaries of a 2D box with x
horizontal and y vertical, boundary conditions of pressure and
velocity which correspond to an ideal-fluid travelling gravity
wave in moderately deep water.
Within the box, the pressure and velocity are computed from the
momentum and energy equations in the usual way.
It is interesting to compare their values with those
corresponding to the ideal-fluid wave, namely: pdif, udif and vdif
If OBSTCL=T, an obstacle appears on the the sea-bed.
If STRA =T, the stresses and strains in it are computed
ENDDIS
Constants used for connecting p1, u1, v1 with ppot, upot and vpot
at the boundaries.
Declarations:
real(pcow,ucow,vcow,pcoe,ucoe,vcoe)
integer(lstepw) ! duration of whole-field patches
boolean(obstcl)
Settings
STRA = F
obstcl = f
obstcl = t
ucow=1.e10; vcow=1.e10; pcow=1.e10
ucoe=1.e10; vcoe=1.e10; pcoe=1.e5
lstepw=2 ! flow field is fixed to potential for first 2 time steps
nx=50 ! discretization in space
ny=50
nx=5
ny=5
lstep=50 ! and time
lstep=5 ! and time
lstep=100
! surface-wave-defining constants
real(lamda,aa,m,h,sig,g,arg,speed,const,twopi,tim1)
g=9.81 ! gravitational acceleration
h=40.0 ! total depth of water
lamda=120 ! wave-length
lamda=1200 ! wave-length
twopi=2.0*3.1416
m=twopi/lamda
arg=m*h
const=tanh(arg)
sig=(g*m*const)**0.5
speed=sig/m
tlast=lamda/speed
tlast=tlast*lstep/100
tim1=tlast/lstep ! time after first time interval
sig
speed
tlast
tim1
aa=3.0 ! amplitude of wave
aa=0.3 ! amplitude of wave
xulast=lamda ! extent of domain in wave-travel direction
yvlast=0.75*h ! domain extent above ocean floor
steady=f
resfac=0.
resref(1)=0.0
endit(1) =0.0
resref(3)=0.0
endit(3) =0.0
resref(5)=0.0
endit(5) =0.0
isg21=lsweep
#unigrid ! use uniform grid
solve(p1,u1,v1)
rho1=1.000e3
rho1=1.
liter(p1)=1000
store(imb1,pcor,min1)
store(u1sl,v1sl,u1ad,v1ad)
(stored var udif is u1-upot )
(stored var vdif is v1-vpot )
(stored var pdif is p1-ppot )
store(pdcx)
relax(u1,falsdt,1.e6)
relax(v1,falsdt,1.e6)
!!!! patches !!!!!
!!!! for inform-created variables !!!!!
patch(whole,phasem,1,nx,1,ny,1,nz,1,lstep) ! for pot, ppot, upot, vpot
! and initial values
!!!! zero-pressure outer patches
patch(p1x,phasem,1,1,1,ny,1,nz,1,lstep) ! west edge
coval(p1x,p1,fixval,0.0)
patch(pny,phasem,1,nx,ny,ny,1,nz,1,lstep) ! north edge
coval(pny,p1,fixval,0.0)
patch(pnx,phasem,nx,nx,1,ny,1,nz,1,lstep) ! east edge
coval(pnx,p1,fixval,0.0)
!!!! next-inner pressure patches
patch(>ppot2x,phasem,2,2,1,ny-1,1,nz,1,lstep) ! near west edge
coval(>ppot2x,p1,pcoe,1)
patch(>ppotny1,phasem,2,nx-1,ny-1,ny-1,1,nz,1,lstep) ! near north edge
coval(>ppotny1,p1,pcoe,1)
patch(>ppotnx1,phasem,nx-1,nx-1,1,ny-1,1,nz,1,lstep) ! near east edge
coval(>ppotnx1,p1,pcoe,1)
!!!! velocity patches near boundaries
patch(>upot1x,phasem,1,1,1,ny,1,nz,1,lstep) ! near west edge
coval(>upot1x,u1,ucoe,1)
patch(>vpot1x,phasem,1,1,1,ny-1,1,nz,1,lstep)
coval(>vpot1x,v1,vcoe,1)
patch(>upotny,phasem,1,nx,ny,ny,1,nz,1,lstep) ! near north edge
coval(>upotny,u1,ucoe,1)
patch(>vpotny1,phasem,1,nx,ny-1,ny-1,1,nz,1,lstep)
coval(>vpotny1,v1,vcoe,1) ! edges
patch(>upotnx1,phasem,nx-1,nx-1,1,ny-1,1,nz,1,lstep) ! near east edge
coval(>upotnx1,u1,ucoe,1)
patch(>vpotnx,phasem,nx,nx,1,ny-1,1,nz,1,lstep)
coval(>vpotnx,v1,vcoe,1)
!!!! whole-field patches for lstepw time steps
patch(>ppotw,phasem,1,nx,1,ny,1,nz,1,lstepw)
coval(>ppotw,p1,pcow ,1)
patch(>upotw,phasem,1,nx,1,ny,1,nz,1,lstepw)
coval(>upotw,u1,ucow,1) ! whole domain
patch(>vpotw,phasem,1,nx,1,ny-1,1,nz,1,lstepw)
coval(>vpotw,v1,vcow,1) ! whole domain
! defining the velocity potential and derived quantities
char(form) ! declare character variable
Formula for the potential as function of space and time
form=aa*(cosh(m*yg))*cos(m*xg-sig*tim)
(stored var pot is :form:)
form=aa*(cosh(m*yg))*cos(m*xg-sig*tim1) ! note use of tim1
(initial of pot at whole is :form:)
Formula for the potential-derived u velocity =
- d pot/dx Note: d/dx of cos mx = - m*sin ; so plus is OK
form = aa*m*(cosh(m*yg))*sin(m*xu-sig*tim) ! for all times
(stored var upot at whole is :form:)
form = aa*m*(cosh(m*yg))*sin(m*xu-sig*tim1) ! note use of tim1 here
(initial of upot at whole is :form:)
(initial of u1 at whole is :form:)
Formula for the potential-derived v velocity, =
- d pot / dy Note: d/dy of cosh my is sinh my ; so minus is appropriate
form = - aa*m*(sinh(m*yv))*cos(m*xg-sig*tim)
(stored var vpot at whole is :form:)
form = - aa*m*(sinh(m*yv))*cos(m*xg-sig*tim1) ! note use of tim1
(initial of vpot at whole is :form:)
(initial of v1 at whole is :form:)
Formula for potential-derived pressure =
d pot / dt Note: d/dt of cos mx - sig t = -(-sig)*sin ; so plus is OK
const=sig*aa*rho1
form = const*(cosh(m*yg))*sin(m*xg-sig*tim)
(stored var ppot at whole is :form:)
form = const*(cosh(m*yg))*sin(m*xg-sig*tim1) ! note use of tim1
(initial of ppot at whole is :form:)
(initial of p1 at whole is :form:)
if(nz.eq.1) then
idispa=lstep/50 ! for dumping to parphi
endif
text(p :pcow: :pcoe: u :ucow: :ucoe:
obstacle to u1 on floor
patch(obstcl1,phasem,nx/2,nx/2,1,ny/2,1,1,1,lstep)
patch(obstcl1,phasem,nx/2,nx/2,1,ny/2,1,1,1,lstep)
coval(obstcl1,u1,1.e10,0.0)
structure on ocean floor
__________
C ! ____ !
B ! !//! !
! !//! !
! !//! !
! !//! !
A ! !//! !
!//! !
--------------------------!//!--!-----------------------------------
1 2 3 4
integer(ix1,ix2,ix3,ix4,iyA,iyB,iyC)
ix1=nx/2
ix2=ix1+2
ix3=ix2
ix4=ix3+2
iyA=5
iyB=iyA+10
iyC=iyB+5
store(prps)
fiinit(prps)= -1
if(obstcl) then
iya=1
patch(wall11,inival,ix1,ix1,iyA,iyC,1,1,1,lstep)
init(wall11,prps,0,100)
patch(wall23,inival,ix2,ix3,1,iyB,1,1,1,lstep)
init(wall23,prps,0,100)
patch(wall44,inival,ix4,ix4,1,iyC,1,1,1,lstep)
init(wall44,prps,0,100)
patch(topC,inival,ix1,ix4,iyC+1,iyC+1,1,1,1,lstep)
init(topC,prps,0,100)
patch(bottom,phasem,ix3,ix4,1,1,1,1,1,lstep)
coval(bottom,p1,fixval,0.0)
endif
#maxabs
# the number of spedat special-print commands
SPEDAT(PRINT,NUMBER,I,3)
# the command to print minimum and maximum values
SPEDAT(PRINT,COMMAND1,C,MINMAX_p1)
SPEDAT(PRINT,COMMAND2,C,MINMAX_u1)
SPEDAT(PRINT,COMMAND3,C,MINMAX_v1)
SPEDAT(PRINT,COMMAND4,C,MINMAX_v1sl)
SPEDAT(PRINT,COMMAND5,C,MINMAX_v1ad)
SPEDAT(PRINT,COMMAND6,C,MINMAX_imb1)
SPEDAT(PRINT,COMMAND7,C,MINMAX_pcor)
-----------------------------------------
if(stra) then
Additions for stress
STORE(EPSZ,EPSY,EPSX,STRZ,STRY,STRX,DRH1,DVO1,ENUL,DEN1,PRPS)
PATCH(BOTTOM11,SOUTH , ix1, ix1, 1, 1, 1, 1, 1, 50)
COVAL(BOTTOM11,U1 , 1.000000E+10, 0.000000E+00)
COVAL(BOTTOM11,V1 , 1.000000E+10, 0.000000E+00)
PATCH(BOTTOM23,SOUTH , ix2, ix3, 1, 1, 1, 1, 1, 50)
COVAL(BOTTOM23,U1 , 1.000000E+10, 0.000000E+00)
COVAL(BOTTOM23,V1 , 1.000000E+10, 0.000000E+00)
PATCH(BOTTOM44,SOUTH , ix4, ix4, 1, 1, 1, 1, 1, 50)
COVAL(BOTTOM44,U1 , 1.000000E+10, 0.000000E+00)
COVAL(BOTTOM44,V1 , 1.000000E+10, 0.000000E+00)
SPEDAT(SET,RLXFAC,RLXU1D,R,7.50000E-01)
SPEDAT(SET,LONGNAME,P1,C,PRESSURE_OR_DILATATION)
SPEDAT(SET,LONGNAME,ENUL,C,VISCOSITY_OR_LPLUS2G*1.E11)
SPEDAT(SET,LONGNAME,DRH1,C,COMPRESSIBILITY_OR_LPLUSG*1.E11)
SPEDAT(SET,LONGNAME,DVO1,C,TH.EXP.CO_R_OR_H*1.E11)
SPEDAT(SET,LONGNAME,U1,C,X-DIRECTION_VELOCITY_OR_DISPLACEM$)
SPEDAT(SET,LONGNAME,U1,C,ENT)
SPEDAT(SET,LONGNAME,EPSX,C,X-DIRECTION_STRAIN)
SPEDAT(SET,LONGNAME,STRX,C,X-DIRECTION_STRESS)
SPEDAT(SET,LONGNAME,V1,C,Y-DIRECTION_VELOCITY_OR_DISPLACEM$)
SPEDAT(SET,LONGNAME,V1,C,ENT)
SPEDAT(SET,LONGNAME,EPSY,C,Y-DIRECTION_STRAIN)
SPEDAT(SET,LONGNAME,STRY,C,Y-DIRECTION_STRESS)
SPEDAT(SET,LONGNAME,U1,C,X-DIRECTION_VELOCITY_OR_DISPLACEM$)
SPEDAT(SET,LONGNAME,U1,C,ENT)
SPEDAT(SET,LONGNAME,V1,C,Y-DIRECTION_VELOCITY_OR_DISPLACEM$)
SPEDAT(SET,LONGNAME,V1,C,ENT)
SPEDAT(SET,LONGNAME,EPSZ,C,Z-DIRECTION_STRAIN)
SPEDAT(SET,LONGNAME,STRZ,C,Z-DIRECTION_STRESS)
SPEDAT(SET,LONGNAME,EPST,C,LINEARTHERMALEXPANSION)
SPEDAT(SET,PRINT,NUMBER,I,7)
SPEDAT(SET,PRINT,COMMAND1,C,MINMAX_P1)
SPEDAT(SET,PRINT,COMMAND2,C,MINMAX_U1)
SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_V1)
SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_strx)
SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_stry)
SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_epsx)
SPEDAT(SET,PRINT,COMMAND3,C,MINMAX_epsy)
endif ! end of stress additions
difcut=0
ixmon=nx/2;iymon=ny/2;izmon=nz/2
tstswp=-1
nxprin=1
nyprin=1
iyprf=1;iyprl=iyc
ixprf=ix1-2;ixprl=ix4-2
ntprin=lstep/5
nprint=1
informbegin
debug t
initial f
store t
formula t
informend
debug=t
lsweep=25
lsweep=10
lsweep=1
lstep=3
ntprin=1
istprf=3
istdb1=3
istdb2=3
iyprf=1
iyprl=2
ixprf=1;ixprl=4
dbindx=t
search=t;idbf0=108591
dbcomp=t
dbsoda=t
dbcmpn=t
dbcmpe=t
dbsol2=t
dbgphi(1)=t
dbgphi(3)=t
dbgphi(5)=t
dbgphi(9)=t
dbflux=t
dbconv=t
stop