RiverVortexTest

From Gerris

(Difference between revisions)
Jump to: navigation, search
Revision as of 01:21, 18 May 2011
Gjrickard (Talk | contribs)

← Previous diff
Revision as of 01:31, 18 May 2011
Gjrickard (Talk | contribs)
(River Tidal Turbine 2D Problem)
Next diff →
Line 41: Line 41:
<source lang=gfs> <source lang=gfs>
- ## Domain of unit size in x and y+# run me using: gerris2D -m
-1 0 GfsSimulation GfsBox GfsGEdge {y = 0.0 x = -0.5} {+
- ## Uniform refinement across the domain+ Define RXMID1 (1.0)
- Refine 7+ Define RYMID1 (0.0)
 + Define rmag1(rx,ry) ( sqrt((rx-RXMID1)*(rx-RXMID1)+(ry-RYMID1)*(ry-RYMID1)) )
 + Define RBIT1V(rx,ry) (rmag1(rx,ry) > 0.0 && rmag1(rx,ry) < 0.15 ? 1 : 0.0)
 + Define RBIT2V(rx,ry) ( (fabs(ry) < 0.4 && rx > -0.4 && rx < 2.4) ? 1 : 0.0)
 + Define RBIT1S(rx,ry) ( exp(-rmag1(rx,ry)/0.01) )
- Time {end = 5.0 dtmax=0.001}+ Define RBIT0V(rx,ry) (rmag1(rx,ry) > 0.0 && rmag1(rx,ry) < 0.35 ? 1 : 0.0)
- ## Initialise single tracer, no advection scheme+ Define LEVMAX 8
- VariableTracer {} R {scheme = none}+ Define LEVMID 8
 + Define LEVMIN 5
 + Define LEVBIG 9
 + Define LEVBIGGER 11
- ## Constant diffusion coefficient+ Define R2 (RBIT2V(rx,ry))
- SourceDiffusion {} R 0.02+ Define R3 (RBIT1V(rx,ry))
 + Define LEVEL1 ( R3 ? LEVMAX : LEVMIN )
 + Define LEVEL2 ( R2 ? LEVMAX+1 : LEVEL1)
- ## Initialise tracer to a constant value+# Turbine parameters, diameter is in y direction, width in x direction,
- Init {} {R = 34.45 }+Define DIAMETER 20.
 +Define WIDTH 10.
 +Define Cdbig 12.
 + 
 +#Spacing distance between turbines
 +Define XSPACING DIAMETER*5.
 +Define YSPACING DIAMETER*1.5
 + 
 + 
 +# M2 tidal frequency. The period is 12h25 (44700 seconds) and M2 tidal elevation.
 + Define M2F (2.*M_PI/44700.)
 + Define M2(t) (0.0*cos (M2F*t)+0.01*sin (M2F*t))
 +# Variable U is velocity*depth here, hence set
 +# M2U to be velocity (a number) by depth (P)...
 +#
 + Define M2U(t) (1.0*(P)*sin (M2F*t))
 +# Define tau (-0.8e-4)
 + Define tau (0.0)
 + Define cd (3.e-3 + TT(x,y)*Cdbig)
 + Define theta (1.0)
 + Define onemth (1.0 - theta)
 + Define abfac (dt*cd*Velocity/P)
 + Define bffac (tau*dt)
 + Define onepa (1.+abfac)
 + Define bth (bffac*theta)
 + Define bonemt (bffac*onemth)
 + 
 +# Use the "GfsRiver" model
 +3 2 GfsRiver GfsBox GfsGEdge {} {
 + 
 + Global {
 + /* Dimensional initial potential temperature profile */
 + #define TT(x,y) (fabs(x) < WIDTH/2. && fabs(y) < DIAMETER/2. ? 1. : 0.0)
 + }
 + 
 + # Box size (metres)
 + PhysicalParams { L = 250.0 g = 9.81 }
 + 
 + # I want the centre of the 2nd box at x = 0 so translate by x = 1000
 + GfsMapTransform { tx = -250 ty = 0 tz = 0 }
 + 
 + # End the simulation
 + Time { end = 8640 dtmax = 10 }
 + 
 + # Use terrain module
 + GModule terrain
 + 
 + Refine LEVEL1
 + 
 + VariableTerrain Zb {
 + basename = bathy50Wide
 + path = /data01/marine-physics/rickard/FROMTIM/COSB:/data01/marine-physics/rickard/FROMTIM
 + } { reconstruct = 1 }
 + Init {} {
 + P = MAX (0., -Zb)
 + }
 + 
 + ## remember the initial grid Level...
 + Init{} {MYMINL = Level
 + RSEE1V = RBIT1V(rx,ry)
 + RSEE2V = RBIT2V(rx,ry)
 + CTT = TT(x,y)
 + }
 + 
 + #AdvectionParams { cfl = 1.0 }
 + 
 + AdvectionParams { cfl = 0.5
 + # gradient = gfs_center_sweby_gradient
 + # gradient = gfs_center_superbee_gradient
 + # gradient = gfs_center_van_leer_gradient
 + }
- ApproxProjectionParams { tolerance = 1e-6 } 
ProjectionParams { tolerance = 1e-6 } ProjectionParams { tolerance = 1e-6 }
 + ApproxProjectionParams { tolerance = 1e-6 }
- OutputTime { istep = 10 } stdout+ AdaptGradient { start=360 istep = 1 } {
- OutputTiming { start = end } stdout+ cmax = 0.06
- OutputProjectionStats { istep = 10 } stdout+ cfactor = 2
- OutputSimulation { step = 1.0 } outice_test2d_1-%9.7f.gfs+ maxlevel = LEVEL2
- OutputSimulation { start = end } resice_test2d_1.gfs {}+ minlevel = MYMINL
 + } (P < 1e-4 ? 0. : P + Zb)
- ## Output at a number of points to compare with 
- ## analytic solution 
- GfsOutputLocation { istep = 10 } c1_1dhf_01 { -0.999 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_02 { -0.949 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_03 { -0.899 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_04 { -0.849 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_05 { -0.799 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_06 { -0.749 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_07 { -0.699 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_08 { -0.649 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_09 { -0.599 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_10 { -0.549 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_11 { -0.499 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_12 { -0.449 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_13 { -0.399 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_14 { -0.349 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_15 { -0.299 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_16 { -0.249 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_17 { -0.199 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_18 { -0.149 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_19 { -0.099 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_20 { -0.049 0.0 -0.5 } 
- GfsOutputLocation { istep = 10 } c1_1dhf_21 { -0.001 0.0 -0.5 } 
-}+ # Bottom friction parameterisation semi-implict, plus
 + # plus Coriolis terms added with implicitness parameter theta...
 + #
 + Init { start=360 istep = 1 } {
 + U = ( P > 1e-6 ? (onepa*(U-bonemt*V)-bth*(V+bonemt*U))/(onepa*onepa + bth*bth) : 0.)
 + V = ( P > 1e-6 ? (onepa*(V+bonemt*U)+bth*(U-bonemt*V))/(onepa*onepa + bth*bth) : 0.)
 + }
 + # Try to construct a Vorticity by dividing fluxes by local
 + # depth and applying gradient operators...here in 2D of course...
 + #
 + Init { istep = 1 } {
 + RU = ( P < 1e-4 ? 0.0 : (U/P))
 + RV = ( P < 1e-4 ? 0.0 : (V/P))
 + VORT = (P < 1e-4 ? 0. : (dx("RV") - dy("RU")))
 + }
 + # Vorticity criterion
 + # Takes 1 m/s as normalising maximum velocity
 + AdaptFunction { start = 360 istep = 1 } {
 + cmax = 5e-2
 + cfactor= 2
 + maxlevel = LEVEL2
 + minlevel = MYMINL
 + } (fabs(VORT*dL/1.0))
 +
 + EventBalance { istep = 10 } 0.1
 +
 + OutputTime { istep = 100 } stdout
 +
 + OutputSimulation { step = 100 } boi_wiki_3-%g.gfs
 +
 + OutputSimulation { start = end } boi_wiki_3_end.gfs
 +
 +}GfsBox {
 + left = Boundary {
 + BcDirichlet U M2U(t)
 + BcNeumann P 0
 + }
 + top = Boundary
 + bottom = Boundary
 +}
 +GfsBox { top = Boundary
 + bottom = Boundary
 +}
GfsBox { GfsBox {
- # Gradient condition value is (Flux/Diff Coeff).+ right = Boundary {
- #+ BcDirichlet U M2U(t)
- # If positive value, gradient increases along outward normal+ BcNeumann P 0
- # and implied heat flux is inwards counter to outward normal.+ }
- # If negative value, gradient decreases along outward normal+ top = Boundary
- # and implied heat flux is outwards along outward normal.+ bottom = Boundary
- #+
- left = Boundary { BcNeumann R (0.0063/0.02) }+
- right = Boundary { BcNeumann R (0.0) }+
} }
 +1 2 right
 +2 3 right
 +
</source> </source>

Revision as of 01:31, 18 May 2011

Testing the River version of Gerris. In particular interested in the possibility of tidal turbine like simulations in support of Tim Divett's thesis.

The basic issue is that we note in both the Linear Free Surface and Nonlinear Free Surface versions of Gerris that vortices generated from wakes caused by localised regions of high bottom friction (mimicing the presence of a turbine, say) can sometimes become unstable, i.e., the vortex continues to spin up as it moves downstream leading to (presumably) unphysical magnitudes of flow and fluid height. In the linear free surface version, such vortices persist and can be seen moving with the background flow. In the nonlinear free surface version, these vortices become tighter and more intense, and eventually grind the simulation to a halt via excessive magnitudes of flow and amplitude.

This issue is also relevant to a whole series of linear free surface tidal simulations around headlands and bays around New Zealand. In these simulations an increase in spatial resolution reveals the generation of many eddies that appear to be shed in a physical way via interactions with the coastline. As above, these eddies persist in the background tidal flow. They do not appear to be unstable (ie the runs keep going), but the question remains about the fidelity and actuality of these eddies (and to test this we need more field observations, and perhaps some more thought about the model).

Both models have been tested against the cases present by Signell and Geyer (1991) and both seem to do a good job of reproducing many of the features seen. Thus the underlying schemes seem to be correct, at least for the scales of interest in these tests.

Perhaps it's not surprising that these features arise. After all the models in question here are solving shallow water variants of reality, so that the likely possibility of 3D processes in the water column influencing the vortex behaviour is absent. For the shallow water the only explicit damping is via a bottom friction term, and the only implicit damping is through truncation terms in the advection scheme; these may be missing important physics that allow for a more complete representation of the vortex behaviour.

So, here show some examples of the unstable vortex behaviour in the River code. This can perhaps lead to a better understanding of the limitations of the shallow water models, and how far we should be pushing them to look at the (presumably) complex tidal interactions with coastline and in-situ structures.

River Tidal Turbine 2D Problem

An analytic 1D heat flux solution compared with a Gerris2D run. Run out to t=5, taking about half an hour of single processor time.

# run me using: gerris2D -m
 
Define RXMID1 (1.0)
Define RYMID1 (0.0)
Define rmag1(rx,ry) ( sqrt((rx-RXMID1)*(rx-RXMID1)+(ry-RYMID1)*(ry-RYMID1)) )
Define RBIT1V(rx,ry) (rmag1(rx,ry) > 0.0 && rmag1(rx,ry) < 0.15 ? 1 : 0.0)
Define RBIT2V(rx,ry) ( (fabs(ry) < 0.4 && rx > -0.4 && rx < 2.4) ? 1 : 0.0)
Define RBIT1S(rx,ry) ( exp(-rmag1(rx,ry)/0.01) )
 
Define RBIT0V(rx,ry) (rmag1(rx,ry) > 0.0 && rmag1(rx,ry) < 0.35 ? 1 : 0.0)
 
Define LEVMAX 8
Define LEVMID 8
Define LEVMIN 5
Define LEVBIG 9
Define LEVBIGGER 11
 
Define R2 (RBIT2V(rx,ry))
Define R3 (RBIT1V(rx,ry))
Define LEVEL1 ( R3 ? LEVMAX : LEVMIN )
Define LEVEL2 ( R2 ? LEVMAX+1 : LEVEL1)
 
# Turbine parameters, diameter is in y direction, width in x direction,
Define DIAMETER 20.
Define WIDTH 10.
Define Cdbig 12.
 
#Spacing distance between turbines
Define XSPACING DIAMETER*5.
Define YSPACING DIAMETER*1.5
 
 
# M2 tidal frequency. The period is 12h25 (44700 seconds) and M2 tidal elevation.
Define M2F (2.*M_PI/44700.)
Define M2(t) (0.0*cos (M2F*t)+0.01*sin (M2F*t))
# Variable U is velocity*depth here, hence set
# M2U to be velocity (a number) by depth (P)...
#
Define M2U(t) (1.0*(P)*sin (M2F*t))
# Define tau (-0.8e-4)
Define tau (0.0)
Define cd (3.e-3 + TT(x,y)*Cdbig)
Define theta (1.0)
Define onemth (1.0 - theta)
Define abfac (dt*cd*Velocity/P)
Define bffac (tau*dt)
Define onepa (1.+abfac)
Define bth (bffac*theta)
Define bonemt (bffac*onemth)
 
# Use the "GfsRiver" model
3 2 GfsRiver GfsBox GfsGEdge {} {
 
Global {
/* Dimensional initial potential temperature profile */
#define TT(x,y) (fabs(x) < WIDTH/2. && fabs(y) < DIAMETER/2. ? 1. : 0.0)
}
 
# Box size (metres)
PhysicalParams { L = 250.0 g = 9.81 }
 
# I want the centre of the 2nd box at x = 0 so translate by x = 1000
GfsMapTransform { tx = -250 ty = 0 tz = 0 }
 
# End the simulation
Time { end = 8640 dtmax = 10 }
 
# Use terrain module
GModule terrain
 
Refine LEVEL1
 
VariableTerrain Zb {
basename = bathy50Wide
path = /data01/marine-physics/rickard/FROMTIM/COSB:/data01/marine-physics/rickard/FROMTIM
} { reconstruct = 1 }
Init {} {
P = MAX (0., -Zb)
}
 
## remember the initial grid Level...
Init{} {MYMINL = Level
RSEE1V = RBIT1V(rx,ry)
RSEE2V = RBIT2V(rx,ry)
CTT = TT(x,y)
}
 
#AdvectionParams { cfl = 1.0 }
 
AdvectionParams { cfl = 0.5
# gradient = gfs_center_sweby_gradient
# gradient = gfs_center_superbee_gradient
# gradient = gfs_center_van_leer_gradient
}
 
ProjectionParams { tolerance = 1e-6 }
ApproxProjectionParams { tolerance = 1e-6 }
 
AdaptGradient { start=360 istep = 1 } {
cmax = 0.06
cfactor = 2
maxlevel = LEVEL2
minlevel = MYMINL
} (P < 1e-4 ? 0. : P + Zb)
 
 
# Bottom friction parameterisation semi-implict, plus
# plus Coriolis terms added with implicitness parameter theta...
#
Init { start=360 istep = 1 } {
U = ( P > 1e-6 ? (onepa*(U-bonemt*V)-bth*(V+bonemt*U))/(onepa*onepa + bth*bth) : 0.)
V = ( P > 1e-6 ? (onepa*(V+bonemt*U)+bth*(U-bonemt*V))/(onepa*onepa + bth*bth) : 0.)
}
 
# Try to construct a Vorticity by dividing fluxes by local
# depth and applying gradient operators...here in 2D of course...
#
Init { istep = 1 } {
RU = ( P < 1e-4 ? 0.0 : (U/P))
RV = ( P < 1e-4 ? 0.0 : (V/P))
VORT = (P < 1e-4 ? 0. : (dx("RV") - dy("RU")))
}
# Vorticity criterion
# Takes 1 m/s as normalising maximum velocity
AdaptFunction { start = 360 istep = 1 } {
cmax = 5e-2
cfactor= 2
maxlevel = LEVEL2
minlevel = MYMINL
} (fabs(VORT*dL/1.0))
 
EventBalance { istep = 10 } 0.1
 
OutputTime { istep = 100 } stdout
 
OutputSimulation { step = 100 } boi_wiki_3-%g.gfs
 
OutputSimulation { start = end } boi_wiki_3_end.gfs
 
}GfsBox {
left = Boundary {
BcDirichlet U M2U(t)
BcNeumann P 0
}
top = Boundary
bottom = Boundary
}
GfsBox { top = Boundary
bottom = Boundary
}
GfsBox {
right = Boundary {
BcDirichlet U M2U(t)
BcNeumann P 0
}
top = Boundary
bottom = Boundary
}
1 2 right
2 3 right

Top frame compares space-time solution from Gerris2D (the red dashes) and from the analytic Laplace transform solution.

Green crosses indicate x-positions of sampled output.

Lower frame plots changes in concentration as a function of time at the indicated x-positions. The dashed lines are from Gerris2D on top of the Laplace Transform solutions.

Comparison of 1D Heat Flux solutions.
Gerris
Enlarge
Gerris
Personal tools
communication