# RiverVortexTest

### From Gerris

Revision as of 09:41, 19 May 2011Popinet (Talk | contribs) (→River Tidal Turbine 2D Sample Output) ← Previous diff |
Revision as of 12:13, 19 May 2011Popinet (Talk | contribs) (Added 2D NS results) Next diff → |
||

Line 258: |
Line 258: | ||

vorticity criterion which certainly seems to work as expected. Thus the | vorticity criterion which certainly seems to work as expected. Thus the | ||

evolving vortices are resolved to the finest scale set. | evolving vortices are resolved to the finest scale set. | ||

+ | |||

+ | == 2D incompressible Euler results == | ||

+ | |||

+ | Given that the Froude number for this example is only Fr = 1./sqrt(g*H) = 0.045, it is entirely justified to use a rigid-lid approximation. The problem thus becomes depth-independent an can be solved using the incompressible Euler (no viscous terms) solver ([[GfsSimulation]]) with the same "bottom friction" parameterisation. | ||

+ | |||

+ | The parameter file is very similar to the GfsRiver one above | ||

+ | |||

+ | <source gfs> | ||

+ | # 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) (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/50.) | ||

+ | Define bffac (tau*dt) | ||

+ | Define onepa (1.+abfac) | ||

+ | Define bth (bffac*theta) | ||

+ | Define bonemt (bffac*onemth) | ||

+ | |||

+ | # Use the "GfsRiver" model | ||

+ | 3 2 GfsSimulation 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 } | ||

+ | |||

+ | # 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 } | ||

+ | |||

+ | Refine LEVEL1 | ||

+ | |||

+ | ## remember the initial grid Level... | ||

+ | Init{} {MYMINL = Level | ||

+ | RSEE1V = RBIT1V(rx,ry) | ||

+ | RSEE2V = RBIT2V(rx,ry) | ||

+ | CTT = TT(x,y) | ||

+ | } | ||

+ | |||

+ | # Bottom friction parameterisation semi-implict, plus | ||

+ | # plus Coriolis terms added with implicitness parameter theta... | ||

+ | # | ||

+ | Init { start=360 istep = 1 } { | ||

+ | U = (onepa*(U-bonemt*V)-bth*(V+bonemt*U))/(onepa*onepa + bth*bth) | ||

+ | V = (onepa*(V+bonemt*U)+bth*(U-bonemt*V))/(onepa*onepa + bth*bth) | ||

+ | } | ||

+ | |||

+ | # Try to construct a Vorticity by dividing fluxes by local | ||

+ | # depth and applying gradient operators...here in 2D of course... | ||

+ | # | ||

+ | Init { istep = 1 } { | ||

+ | VORT = (dx("V") - dy("U")) | ||

+ | } | ||

+ | # 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 = 10 } stderr | ||

+ | OutputProjectionStats { istep = 10 } stderr | ||

+ | |||

+ | # OutputSimulation { istep = 10 } stdout | ||

+ | |||

+ | OutputSimulation { step = 100 } boi_wiki_3-%g.gfs | ||

+ | |||

+ | OutputSimulation { start = end } boi_wiki_3_end.gfs | ||

+ | |||

+ | GModule gfsview | ||

+ | OutputView { step = 20 } { ppm2mpeg > vort.mpg } vort.gfv | ||

+ | } | ||

+ | GfsBox { | ||

+ | left = Boundary { | ||

+ | BcDirichlet U M2U(t) | ||

+ | BcNeumann P 0 | ||

+ | } | ||

+ | top = Boundary | ||

+ | bottom = Boundary | ||

+ | } | ||

+ | GfsBox { top = Boundary | ||

+ | bottom = Boundary | ||

+ | } | ||

+ | GfsBox { | ||

+ | right = BoundaryOutflow | ||

+ | top = Boundary | ||

+ | bottom = Boundary | ||

+ | } | ||

+ | 1 2 right | ||

+ | 2 3 right | ||

+ | </source> | ||

+ | |||

+ | The results at t = 3500 look like this. Of course, these results will be grid-dependent since no explicit viscosity has been added. | ||

+ | |||

+ | {| align="center" | ||

+ | |+ Zonal speed (colour), height (contours). | ||

+ | |- | ||

+ | | [[Image:Euler_3500_0.png|thumb|400px|center|2D Euler]] | ||

+ | | [[Image:Sweby_3500_0.png|thumb|400px|center|Saint-Venant Sweby]] | ||

+ | |} | ||

+ | |||

+ | {| align="center" | ||

+ | |+ Height (colour and contours). | ||

+ | |- | ||

+ | | [[Image:Euler_3500_1.png|thumb|400px|center|2D Euler]] | ||

+ | | [[Image:Sweby_3500_1.png|thumb|400px|center|Saint-Venant Sweby]] | ||

+ | |} | ||

+ | |||

+ | {| align="center" | ||

+ | |+ Vorticity (colour) and height (contours). | ||

+ | |- | ||

+ | | [[Image:Euler_3500_2.png|thumb|400px|center|2D Euler]] | ||

+ | | [[Image:Sweby_3500_2.png|thumb|400px|center|Saint-Venant Sweby]] | ||

+ | |} | ||

+ | |||

== Ocean Simulation Output == | == Ocean Simulation Output == |

## Revision as of 12:13, 19 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.

## Contents |

## River Tidal Turbine 2D Problem

A typical 2D run set up with constant bathymetry of 50 metres depth, with an increased region of bottom friction defined by TT(x,y) in the centre of the domain. M2U(t) drives the tidal flow at the boundary, anticipating a maximum input speed of 1 m/s.

The change of relevance here is the choice of advection scheme. The default (Minmod) produces a stable simulation, with a smooth wake behind the turbine. Change that to Sweby, and the wake becomes unstable, leading to vortices that then continue to deepen and amplify.

The question is which is "correct"?

Run out to around t=5000, taking about the same number of cpu seconds on 6 processor elements.

# 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

## River Tidal Turbine 2D Sample Output

Solutions at 3500 seconds into each run using the default (left column) and Sweby (right column) advection schemes.

As perhaps expected the default (minmod) advection gives a smoother wake, with no hint at this time of wake instability and vortex generation.

In contrast with Sweby advection instability is clear and looks (to begin with) physical. However, as these frames show, there is a fairly 'singular' vortex formed, with deepening amplitude and associated flow speeds. This then later generates a whole raft of large amplitude waves over the whole domain, killing the simulation; this process is clearly unstable.

Note that adaptivity is turned on. Initially the grid is refined to the finest level in a zone around the turbine, with relatively coarse cells elsewhere. This region is then allowed to follow the evolution of the wakes through a vorticity criterion which certainly seems to work as expected. Thus the evolving vortices are resolved to the finest scale set.

## 2D incompressible Euler results

Given that the Froude number for this example is only Fr = 1./sqrt(g*H) = 0.045, it is entirely justified to use a rigid-lid approximation. The problem thus becomes depth-independent an can be solved using the incompressible Euler (no viscous terms) solver (GfsSimulation) with the same "bottom friction" parameterisation.

The parameter file is very similar to the GfsRiver one above

# 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) (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/50.)

Define bffac (tau*dt)

Define onepa (1.+abfac)

Define bth (bffac*theta)

Define bonemt (bffac*onemth)

# Use the "GfsRiver" model

3 2 GfsSimulation 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 }

# 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 }

Refine LEVEL1

## remember the initial grid Level...

Init{} {MYMINL = Level

RSEE1V = RBIT1V(rx,ry)

RSEE2V = RBIT2V(rx,ry)

CTT = TT(x,y)

}

# Bottom friction parameterisation semi-implict, plus

# plus Coriolis terms added with implicitness parameter theta...

#

Init { start=360 istep = 1 } {

U = (onepa*(U-bonemt*V)-bth*(V+bonemt*U))/(onepa*onepa + bth*bth)

V = (onepa*(V+bonemt*U)+bth*(U-bonemt*V))/(onepa*onepa + bth*bth)

}

# Try to construct a Vorticity by dividing fluxes by local

# depth and applying gradient operators...here in 2D of course...

#

Init { istep = 1 } {

VORT = (dx("V") - dy("U"))

}

# 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 = 10 } stderr

OutputProjectionStats { istep = 10 } stderr

# OutputSimulation { istep = 10 } stdout

OutputSimulation { step = 100 } boi_wiki_3-%g.gfs

OutputSimulation { start = end } boi_wiki_3_end.gfs

GModule gfsview

OutputView { step = 20 } { ppm2mpeg > vort.mpg } vort.gfv

}

GfsBox {

left = Boundary {

BcDirichlet U M2U(t)

BcNeumann P 0

}

top = Boundary

bottom = Boundary

}

GfsBox { top = Boundary

bottom = Boundary

}

GfsBox {

right = BoundaryOutflow

top = Boundary

bottom = Boundary

}

1 2 right

2 3 right

The results at t = 3500 look like this. Of course, these results will be grid-dependent since no explicit viscosity has been added.

## Ocean Simulation Output

Using the GfsOcean simulation we get some similar super vortices generated as the tide reverses back through an array of turbines. The super vortex (red area above the array in figure below) forms to the right of the array at the reversing of the tide and passes near the top right turbine in the array which whips across to its current position. It grows in magnitude as it moves. I don't see this in a simulation of a single turbine.

## Summary

This form of instability is not peculiar to the nonlinear free surface code, but the nonlinearity feeds the effect of the vortex to the wider field, causing a model crash. The growth of the vortex is similar in the linear free surface code (might need an example from Tim) but it doesn't generally lead to a model crash.

I suspect that we're expecting too much from these versions of the model, ie, shallow water with their own advection scheme variants. It may well be that a full resolution will have to be in 3D (and indeed Sebastien has tried that and found how tricky and expensive that is!).

Underlying the problem remains the fundamental issue; what is the eddy generation potential of tidal flow interacting with coastline? The Gerris models suggest that there are potentially a lot of eddies out there to be found; however the fact that the model response is ultimately unstable suggests that these models might be over-estimating that potential.