-
Notifications
You must be signed in to change notification settings - Fork 0
/
relaxedGreenNaghdi.jl
291 lines (253 loc) · 10.4 KB
/
relaxedGreenNaghdi.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
export relaxedGreenNaghdi
"""
relaxedGreenNaghdi(param;kwargs)
Define an object of type `AbstractModel` in view of solving the initial-value problem for
a relaxed Green-Naghdi model proposed by [N. Favrie and S. Gavrilyuk](https://doi.org/10.1088/1361-6544/aa712d)
or [C. Escalante, M. Dumbser and M. Castro](https://doi.org/10.1016/j.jcp.2019.05.035)
and [G. Richard](https://doi.org/10.1016/j.euromechflu.2021.05.011).
# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- the relaxation parameter `a`;
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.
## Optional keyword arguments
- `FG`: if `true` (default is `false`), compute the Favrie-Gavrilyuk model, otherwise compute the Escalante-Dumbser-Castro model;
- `id`: `∈{0,1,2}` and represent the level of preparation of the initial data (default is `1`);
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `iterative`: solve the elliptic problem (to construct initial data) through GMRES if `true`, LU decomposition if `false` (default is `true`);
- `precond`: use a (left) preconditioner for GMRES if `true` (default), choose `precond` as the preconditioner if provided;
- `gtol`: relative tolerance of the GMRES algorithm (default is `1e-14`);
- `restart`: the corresponding option of the GMRES algorithm (default is `100`);
- `maxiter`: the corresponding option of GMRES (default is `nothing`);
- `ktol`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Favrie-Gavrilyuk"` if `FG==true`, `"Escalante-Dumbser-Castro"` otherwise);
# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `relaxedGreenNaghdi.f!` to be called in explicit time-integration solvers;
2. a function `relaxedGreenNaghdi.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `relaxedGreenNaghdi.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
- `η` is the values of surface deformation at collocation points `x`;
- `v` is the derivative of the trace of the velocity potential at `x`;
4. additionally, a handy function `relaxedGreenNaghdi.mapfrofull` which from data matrix returns the Tuple of real vectors `(η,v,u,p,w)`, where
- `u` corresponds to the layer-averaged horizontal velocity.
- `p` corresponds to the relaxed (artificial) layer-averaged non-hydrostatic pressure;
- `w` corresponds to the relaxed (artificial) layer-averaged vertical velocity.
"""
mutable struct relaxedGreenNaghdi <: AbstractModel
label :: String
f! :: Function
mapto :: Function
mapfro :: Function
mapfrofull :: Function
info :: String
function relaxedGreenNaghdi(param::NamedTuple; FG = false,
id = 1,
mesh = Mesh(param),
dealias = 0,
ktol = 0,
iterate = true,
gtol = 1e-14,
precond = true,
restart = nothing,
maxiter = nothing,
label = nothing
)
# Set up
a = param.a
μ = param.μ
ϵ = param.ϵ
if isnothing(maxiter) maxiter = mesh.N end
if isnothing(restart) restart = min(20,mesh.N) end
if isnothing(label)
if FG == true
label = "Favrie-Gavrilyuk"
else
label = "Escalante-Dumbser-Castro"
end
end
# Print information
info = "$label model.\n"
info *= "├─Relaxation parameter a=$a.\n"
info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ.\n"
if dealias == 0
info *= "├─No dealiasing. "
else
info *= "├─Dealiasing with Orszag's rule adapted to power $(dealias + 1) nonlinearity. "
end
if ktol == 0
info *= "No Krasny filter. "
else
info *= "Krasny filter with tolerance $ktol."
end
if iterate == true
if precond == false out="out" else out="" end
info *= "\n└─Elliptic problem solved with GMRES method with$out preconditioning, \
tolerance $gtol, maximal number of iterations $maxiter, restart after $restart iterations \
(consider `iterate=false` for non-iterative method). "
else
info *= "\n└─Elliptic problem solved with standard LU factorization \
(consider `iterate=true` for faster results). "
end
info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."
# Pre-allocate useful data
k = mesh.k
x = mesh.x
x₀ = mesh.x[1]
∂ₓ = 1im * k
if precond == true
Precond = Diagonal( 1 .+ μ/3*k.^2 )
elseif precond == false
Precond = Diagonal( ones(size(k)) )
else
Precond = precond
end
if dealias == 0
Π⅔ = ones(size(k)) # no dealiasing (Π⅔=Id)
else
K = (mesh.kmax-mesh.kmin)/(2+dealias)
Π⅔ = abs.(k) .<= K # Dealiasing low-pass filter
end
FFT = exp.(-1im*k*(x.-x₀)');
IFFT = exp.(1im*k*(x.-x₀)')/length(x);
Dx = IFFT * Diagonal( Π⅔ .* ∂ₓ)* FFT
Id = Diagonal(ones(size(x)));
h = zeros(Complex{Float64}, mesh.N)
u, p, w, fftη, fftu, fftv, fftp, fftw, hFG = (similar(h),).*ones(9)
L = similar(FFT)
# Evolution equations are ∂t U = f(U)
function f!(U;a=a)
fftη .= U[:,1]
h .= 1 .+ ϵ*ifft(fftη)
fftu .= U[:,2]; u.= ifft(fftu);
fftp .= U[:,3]; p.= ifft(fftp);
fftw .= U[:,4]; w.= ifft(fftw);
if FG == true
hFG .= h.*(1 .- ϵ*p/a)
else
hFG .= h
end
U[:,1] .= -∂ₓ.*Π⅔.*fft(h .* u)
U[:,2] .= -∂ₓ.*Π⅔.*(fftη .+ ϵ/2 * fft( u.^2) ) - a*μ * Π⅔.*fft(1 ./h .* ifft(∂ₓ.* fft( hFG.*p ) ) )
U[:,3] .= -a*Π⅔.*fft((2*w.+hFG.*ifft(∂ₓ.*fftu))./h) - ϵ*Π⅔.*fft( u .* ifft(∂ₓ.* fftp) )
U[:,4] .= 3/2*a*Π⅔.*fft(p./h) - ϵ*Π⅔.*fft( u .* ifft(∂ₓ.* fftw ) )
U[abs.(U).< ktol ].=0
end
# Build raw data from physical data.
# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
function mapto(data::InitialData)
fftη .= fft(data.η(x))
fftv .= fft(data.v(x))
h .= 1 .+ ϵ*ifft(fftη)
if iterate == false
L .= Id - μ/3 * Diagonal(Π⅔) * FFT * Diagonal( 1 ./h ) * Dx * Diagonal( h.^3 ) * IFFT * Diagonal( Π⅔ .* ∂ₓ )
fftu .= L \ fftv
else#if iterate == true
function LL(hatu)
hatu- μ/3 *Π⅔.*fft( 1 ./h .* ifft( Π⅔ .* ∂ₓ .*fft( h.^3 .* ifft( Π⅔ .* ∂ₓ .* hatu ) ) ) )
end
fftu .= gmres( LinearMap(LL, length(h); issymmetric=false, ismutating=false) , fftv ;
restart = restart, maxiter = maxiter, Pl = Precond, reltol = gtol )
end
if id >= 2
if iterate == false
L .= FFT* Diagonal( 3 ./h.^3 ) * IFFT - μ * Diagonal( Π⅔ .* ∂ₓ ) * FFT * Diagonal( 1 ./h ) * IFFT * Diagonal( Π⅔ .* ∂ₓ )
fftp .= fft(ifft(Π⅔ .*(L \ (Π⅔ .* ( ∂ₓ.^2 .*fftη + 2 * fft(ifft(Π⅔ .* ∂ₓ.* fftu).^2)))) / a)./h)
elseif iterate == true
function ll(hatu)
Π⅔.*fft( 3 ./h.^3 .* ifft(hatu))- μ * Π⅔ .* ∂ₓ .*fft( 1 ./h .* ifft( Π⅔ .*∂ₓ .* hatu ) )
end
fftp .= fft(ifft(Π⅔ .* gmres( LinearMap(ll, length(h); issymmetric=false, ismutating=false) , Π⅔ .* ( ∂ₓ.^2 .*fftη + 2 * fft(ifft(Π⅔ .* ∂ₓ.* fftu).^2)) ;
restart = restart, maxiter = maxiter, Pl = Precond, reltol = gtol ) / a) ./h)
elseif iterate == 1/2
U₀ = [Π⅔ .* fftη Π⅔ .* fftu fft(zero(x)) Π⅔ .* (-1/2*fft(h.*ifft(∂ₓ.*fftu)))]
dt = param.dt.^2
# Calcul U(dt)
U0 = copy(U₀)
f!( U0 ; a = a^2)
U1 = copy(U0)
U0 .= U₀ .+ dt/2 .* U1
f!( U0 ; a = a^2)
U2 = copy(U0)
U0 .= U₀ .+ dt/2 .* U2
f!( U0 ; a = a^2)
U3 = copy(U0)
U0 .= U₀ .+ dt .* U3
f!( U0 ; a = a^2)
U4 = copy(U0)
U₊ = U₀ + dt/6 .* (U1 + 2*U2 + 2*U3 + U4 )
# Calcul U(-dt)
U0 = copy(U₀)
f!( U0 ; a = a^2)
U1 = -copy(U0)
U0 .= U₀ .+ dt/2 .* U1
f!( U0 ; a = a^2)
U2 = -copy(U0)
U0 .= U₀ .+ dt/2 .* U2
f!( U0 ; a = a^2)
U3 = -copy(U0)
U0 .= U₀ .+ dt .* U3
f!( U0 ; a = a^2)
U4 = -copy(U0)
U₋ = U₀ + dt/6 .* (U1 + 2*U2 + 2*U3 + U4 )
d2th = ifft(U₊[:,1]+U₋[:,1]-2*U₀[:,1])/dt^2
dth = ifft(U₊[:,1]-U₋[:,1])/(2*dt)
dtdxh = ifft(∂ₓ.*U₊[:,1]-∂ₓ.*U₋[:,1])/(2*dt)
d2xh = ifft(∂ₓ.*∂ₓ.*U₀[:,1])
dtu = ifft(U₊[:,2]-U₋[:,2])/(2*dt)
dxu = ifft(∂ₓ.*U₀[:,2])
dxh = ifft(∂ₓ.*U₀[:,1])
u = ifft(U₀[:,2])
# Need to relaod everything since this has been modified by f!
fftη .= fft(data.η(x))
fftv .= fft(data.v(x))
h .= 1 .+ ϵ*ifft(fftη)
fftu .= fft(u)
fftp = 1/3/a*Π⅔ .*fft( h .* (d2th+dtu.*dxh+2*u.*dtdxh+u.*dxu.*dxh+u.*u.*d2xh) )
end
else
fftp.=fft(zero(x))
end
if id == 0
fftw = fft(zero(x))
else
if FG == true
hFG .= h.*(1 .- ϵ*ifft(fftp))
else
hFG .= h
end
fftw.=-1/2*fft(hFG.*ifft(∂ₓ.*fftu))
end
U = [Π⅔ .* fftη Π⅔ .* fftu Π⅔ .* fftp Π⅔ .* fftw]
U[abs.(U).< ktol ].=0
return U
end
# Reconstruct physical variables from raw data
# Return `(η,v,x)`, where
# - `η` is the surface deformation;
# - `v` is the derivative of the trace of the velocity potential;
# - `x` is the vector of collocation points
function mapfro(U)
fftη .= U[:,1]
h .= 1 .+ ϵ*ifft(fftη)
fftu .= U[:,2]
fftv .= fftu - μ/3 *Π⅔.*fft( 1 ./h .* ifft( Π⅔.* ∂ₓ .* fft( h.^3 .* ifft( Π⅔.* ∂ₓ .* fftu ) ) ) )
real(ifft(fftη)),real(ifft(fftv)),mesh.x
end
# Return `(η,v,u,p,w,x)`, where
# - `η` is the surface deformation;
# - `v` is the derivative of the trace of the velocity potential;
# - `u` corresponds to the layer-averaged horizontal velocity.
# - `p` corresponds to the relaxed (artificial) layer-averaged non-hydrostatic pressure;
# - `w` corresponds to the relaxed (artificial) layer-averaged vertical velocity.
# - `x` is the vector of collocation points
function mapfrofull(U)
fftη .= U[:,1]
h .= 1 .+ ϵ*ifft(fftη)
fftu .= U[:,2]
fftv .= fftu - μ/3 *Π⅔.*fft( 1 ./h .* ifft( Π⅔.* ∂ₓ .* fft( h.^3 .* ifft( Π⅔.* ∂ₓ .* fftu ) ) ) )
real(ifft(fftη)),real(ifft(fftv)),real(ifft(fftu)),real(ifft(U[:,3])),real(ifft(U[:,4])),mesh.x
end
new(label, f!, mapto, mapfro, mapfrofull, info )
end
end