1
! Copyright (C) 2006 Imperial College London and others.
3
! Please see the AUTHORS file in the main source directory for a full list
4
! of copyright holders.
6
! Applied Modelling and Computation Group
7
! Department of Earth Science and Engineering
8
! Imperial College London
10
! amcgsoftware@imperial.ac.uk
12
! This library is free software; you can redistribute it and/or
13
! modify it under the terms of the GNU Lesser General Public
14
! License as published by the Free Software Foundation,
15
! version 2.1 of the License.
17
! This library is distributed in the hope that it will be useful,
18
! but WITHOUT ANY WARRANTY; without even the implied warranty of
19
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20
! Lesser General Public License for more details.
22
! You should have received a copy of the GNU Lesser General Public
23
! License along with this library; if not, write to the Free Software
24
! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
29
! This module implements the 2D shallow water equations within the normal
30
! fluidity binary configured from an ordinary .flml, going through the
31
! usual momentum_equation code path. It is
32
! a different implementation than the one in main/Shallow_Water.F90
33
! that has its own binary and schema.
34
module shallow_water_equations
35
use global_parameters, only: OPTION_PATH_LEN
39
use boundary_conditions
45
! Assemble the shallow water continuity equation by adding the time derivative (\eta^{n+1}-\eta^n)/dt
46
! This is based on assemble_1mat_compressible_projection_cg, but a lot simpler:
47
! The density \rho is here the total water depth (bottom depth+fs elevation). Its time-derivative however
48
! is the same as the free surface time derivative. Pressure is g*\eta, so drhodp is simply 1/g
50
subroutine assemble_shallow_water_projection(state, cmc, rhs, dt, theta_pg, theta_divergence, reassemble_cmc_m)
52
! This only works for single material_phase, so just pass me the first state:
53
type(state_type), intent(inout) :: state
54
! the lhs and rhs to add into:
55
type(csr_matrix), intent(inout) :: cmc
56
type(scalar_field), intent(inout) :: rhs
58
real, intent(in) :: dt
59
real, intent(in) :: theta_pg, theta_divergence
60
logical, intent(in) :: reassemble_cmc_m
62
integer, dimension(:), pointer :: test_nodes
63
type(element_type), pointer :: test_shape
64
real, dimension(:), allocatable :: detwei
65
real, dimension(:), allocatable :: delta_p_at_quad
69
type(vector_field), pointer :: coordinate
70
type(scalar_field), pointer :: pressure, old_pressure
72
ewrite(1,*) 'Entering assemble_shallow_water_projection'
76
! only do all this if we need to make cmc (otherwise we'd be adding repeatedly)
77
if(reassemble_cmc_m) then
78
coordinate=> extract_vector_field(state, "Coordinate")
80
pressure => extract_scalar_field(state, "Pressure")
81
old_pressure => extract_scalar_field(state, "OldPressure")
83
call get_option("/physical_parameters/gravity/magnitude", g)
84
! this is the density in front of the du/dt time-derivative which at the
85
! moment is hard-coded to 1.0
88
allocate(detwei(ele_ngi(pressure, 1)), &
89
delta_p_at_quad(ele_ngi(pressure, 1)))
91
do ele=1, element_count(pressure)
93
test_nodes=>ele_nodes(pressure, ele)
95
test_shape => ele_shape(pressure, ele)
97
delta_p_at_quad = ele_val_at_quad(pressure, ele) - ele_val_at_quad(old_pressure, ele)
99
call transform_to_physical(coordinate, ele, detwei=detwei)
101
! Time derivative: pressure correction \Phi that we are solving for is: \Phi=theta_div*theta_pg*(p^{n+1}-p^*)*dt
102
! Thus time derivative lhs is (\eta^{n+1}-\eta^*)/dt = \Phi/(g*dt*dt/theta_div/theta_pg)
103
call addto(cmc, test_nodes, test_nodes, &
104
shape_shape(test_shape, test_shape, detwei)/(rho0*g*dt*dt*theta_divergence*theta_pg))
105
! Time derivative rhs: -(\eta^*-\eta^n)/dt = (p^*-p^n)/(g*dt)
106
call addto(rhs, test_nodes, &
107
-shape_rhs(test_shape, detwei*delta_p_at_quad)/(rho0*g*dt))
111
deallocate(detwei, delta_p_at_quad)
115
end subroutine assemble_shallow_water_projection
117
subroutine assemble_swe_divergence_matrix_cg(ctp_m, state, ct_rhs)
119
! This only works for single material_phase, so just pass me the first state:
120
type(state_type), intent(inout) :: state
122
! the divergence matrix and rhs to be assembled
123
type(block_csr_matrix), pointer :: ctp_m
124
type(scalar_field), intent(inout), optional :: ct_rhs
127
type(mesh_type), pointer :: test_mesh
129
type(vector_field), pointer :: field
131
integer, dimension(:), pointer :: test_nodes, field_nodes
133
real, dimension(:,:,:), allocatable :: ele_mat, ele_mat_bdy
134
type(element_type), pointer :: field_shape, test_shape
135
real, dimension(:,:,:), allocatable :: dfield_t, dtest_t, dbottom_t
136
real, dimension(:), allocatable :: detwei
137
real, dimension(:,:), allocatable :: normal_bdy
139
real, dimension(:), allocatable :: depth_at_quad
140
real, dimension(:,:), allocatable :: depth_grad_at_quad
142
integer :: ele, sele, dim, xdim, ngi
144
type(vector_field), pointer :: coordinate, velocity
145
type(scalar_field), pointer :: pressure, old_pressure, bottom_depth
146
real :: theta, dt, g, rho0
149
logical :: integrate_by_parts
151
integer, dimension(:,:), allocatable :: field_bc_type
152
type(vector_field) :: field_bc
154
ewrite(1,*) 'In assemble_swe_divergence_matrix_cg'
156
coordinate=> extract_vector_field(state, "Coordinate")
158
pressure => extract_scalar_field(state, "Pressure")
159
old_pressure => extract_scalar_field(state, "OldPressure")
160
bottom_depth => extract_scalar_field(state, "BottomDepth")
162
velocity=>extract_vector_field(state, "Velocity")
164
! note that unlike for compressible we adhere to the option under pressure
165
! (unless DG for which we don't have a choice)
166
integrate_by_parts=have_option(trim(pressure%option_path)// &
167
&"/prognostic/spatial_discretisation/integrate_continuity_by_parts") &
168
&.or. have_option(trim(velocity%option_path)// &
169
&"/prognostic/spatial_discretisation/discontinuous_galerkin")
171
ewrite(2,*) "SWE divergence is integrated by parts: ", integrate_by_parts
173
! note that this is different then compressible, as we don't have
174
! a prognostic density equivalent, just treating it as a nonlinear relaxation here
175
call get_option(trim(velocity%option_path)// &
176
&"/prognostic/temporal_discretisation/relaxation", theta)
177
call get_option("/timestepping/timestep", dt)
178
call get_option("/physical_parameters/gravity/magnitude", g)
179
! this is the density in front of the du/dt time-derivative which at the
180
! moment is hard-coded to 1.0
183
test_mesh => pressure%mesh
184
! the field that ctp_m is multiplied with:
187
if(present(ct_rhs)) call zero(ct_rhs)
189
! Clear memory of arrays being designed
192
ngi = ele_ngi(test_mesh, 1)
193
xdim = coordinate%dim
194
allocate(dtest_t(ele_loc(test_mesh, 1), ngi, xdim), &
195
dfield_t(ele_loc(field, 1), ngi, xdim), &
196
dbottom_t(ele_loc(bottom_depth, 1), ngi, xdim), &
197
ele_mat(field%dim, ele_loc(test_mesh, 1), ele_loc(field, 1)), &
199
depth_at_quad(ngi), &
200
depth_grad_at_quad(xdim, ngi))
202
do ele=1, element_count(test_mesh)
204
test_nodes => ele_nodes(test_mesh, ele)
205
field_nodes => ele_nodes(field, ele)
207
test_shape => ele_shape(test_mesh, ele)
208
field_shape => ele_shape(field, ele)
209
call transform_to_physical(coordinate, ele, test_shape, dshape = dtest_t, detwei=detwei)
211
depth_at_quad = (theta*ele_val_at_quad(pressure, ele) + &
212
(1-theta)*ele_val_at_quad(old_pressure, ele))/(rho0*g) + &
213
ele_val_at_quad(bottom_depth, ele)
215
if(integrate_by_parts) then
217
ele_mat = -dshape_shape(dtest_t, field_shape, detwei*depth_at_quad)
221
! transform the field (velocity) derivatives into physical space
222
call transform_to_physical(coordinate, ele, field_shape, dshape=dfield_t)
223
if (bottom_depth%mesh%shape==field_shape) then
225
else if (bottom_depth%mesh%shape==test_shape) then
228
call transform_to_physical(coordinate, ele, ele_shape(bottom_depth, ele), &
232
assert( test_shape==pressure%mesh%shape )
233
depth_grad_at_quad = (theta*ele_grad_at_quad(pressure, ele, dtest_t) + &
234
(1-theta)*ele_grad_at_quad(old_pressure, ele, dtest_t))/(rho0*g) + &
235
ele_grad_at_quad(bottom_depth, ele, dbottom_t)
237
ele_mat = shape_dshape(test_shape, dfield_t, detwei*depth_at_quad) + &
238
shape_shape_vector(test_shape, field_shape, detwei, depth_grad_at_quad)
242
do dim = 1, field%dim
243
call addto(ctp_m, 1, dim, test_nodes, field_nodes, ele_mat(dim,:,:))
247
deallocate(dtest_t, dfield_t, dbottom_t, &
248
ele_mat, detwei, depth_at_quad, &
251
if(integrate_by_parts) then
253
ngi = face_ngi(field,1)
254
xdim = coordinate%dim
255
allocate(detwei(ngi), &
256
depth_at_quad(ngi), &
257
normal_bdy(xdim, ngi))
258
allocate(ele_mat_bdy(field%dim, face_loc(test_mesh, 1), face_loc(field, 1)))
260
assert(surface_element_count(test_mesh)==surface_element_count(field))
261
allocate(field_bc_type(field%dim, surface_element_count(field)))
262
call get_entire_boundary_condition(field, (/ &
266
field_bc, field_bc_type)
268
do sele = 1, surface_element_count(test_mesh)
270
if(any(field_bc_type(:,sele)==2)&
271
.or.any(field_bc_type(:,sele)==3)) cycle
273
test_shape => face_shape(test_mesh, sele)
274
field_shape => face_shape(field, sele)
276
call transform_facet_to_physical(coordinate, sele, &
280
depth_at_quad = (theta*face_val_at_quad(pressure, sele) + &
281
(1-theta)*face_val_at_quad(old_pressure, sele))/(rho0*g) + &
282
face_val_at_quad(bottom_depth, sele)
284
ele_mat_bdy = shape_shape_vector(test_shape, field_shape, &
285
detwei*depth_at_quad, normal_bdy)
287
do dim = 1, field%dim
288
if((field_bc_type(dim, sele)==1).and.present(ct_rhs)) then
289
call addto(ct_rhs, face_global_nodes(test_mesh, sele), &
290
-matmul(ele_mat_bdy(dim,:,:), &
291
ele_val(field_bc, dim, sele)))
293
call addto(ctp_m, 1, dim, face_global_nodes(test_mesh, sele), &
294
face_global_nodes(field, sele), ele_mat_bdy(dim,:,:))
300
call deallocate(field_bc)
301
deallocate(field_bc_type, depth_at_quad, detwei, ele_mat_bdy)
302
deallocate(normal_bdy)
306
end subroutine assemble_swe_divergence_matrix_cg
308
subroutine shallow_water_equations_check_options
310
character(len=*), dimension(1:4), parameter:: forbidden_bc_types = (/ &
311
"free_surface ", "bulk_formulae", &
312
"drag ", "wind_forcing " /)
313
character(len=OPTION_PATH_LEN):: velocity_option_path, pressure_option_path
317
if (.not. have_option("/material_phase/vector_field::Velocity/prognostic/equation::ShallowWater")) return
319
ewrite(2,*) "Checking shallow water options"
321
call get_option("/geometry/dimension", dim)
323
FLExit("With equation type ShallowWater you need a 2D mesh and configuration")
325
if (have_option("/geometry/spherical_earth")) then
326
FLExit("Equation type ShallowWater is not implemented for spherical_earth")
328
if (have_option("/geometry/ocean_boundaries")) then
329
FLExit("Do not specify /geometry/ocean_boundaries with equation type ShallowWater")
331
if (.not. have_option("/physical_parameters/gravity")) then
332
ewrite(0,*) "Missing option /physical_parameters/gravity " // &
333
& "(you may ignore /physical_parameters/vector_field::GravityDirection)"
334
FLExit("With equation type Shallow water you need to specify gravity")
337
! Options under /material_phase
338
if (option_count("/material_phase")/=1) then
339
FLExit("Equation type ShallowWater only works with a single material_phase")
342
! These don't make sense - so let's just forbid them to avoid confusion
343
if (have_option("/material_phase/equation_of_state") .or. &
344
have_option("/material_phase/scalar_field::Density")) then
345
FLExit("With equation type ShallowWater you're not allowed an equation_of_state or Density field")
350
pressure_option_path = "/material_phase/scalar_field::Pressure/prognostic/"
351
if (.not. have_option(pressure_option_path)) then
352
FLExit("With equation type ShallowWater you need a prognostic Pressure field")
355
if (have_option(trim(pressure_option_path)//"solver/iterative_method::cg")) then
356
ewrite(0,*) "For shallow water equations the pressure matrix is asymmetric"
357
ewrite(0,*) 'Therefore you should not use "cg" as the linear solver'
358
ewrite(0,*) 'Use "gmres" instead'
359
FLExit('Cannot use "cg" as linear solver for Pressure with ShallowWater')
362
if (.not. have_option(trim(pressure_option_path)// &
363
"/spatial_discretisation/continuous_galerkin")) then
364
! it might also work with dg, but someone should test that - definitely won't work with cv
365
FLExit('Equation type ShallowWater only work with a continuous galerkin Pressure')
370
velocity_option_path = "/material_phase/vector_field::Velocity/prognostic/"
371
if (.not. have_option(velocity_option_path)) then
372
FLExit("With equation type ShallowWater you need a prognostic Velocity field")
375
call get_option(trim(velocity_option_path)//"/spatial_discretisation/conservative_advection", beta)
377
ewrite(0,*) "The shallow water equations are implemented in non-conservative form"
378
ewrite(0,*) "To be consistent with that the velocity advection term should be in non-conservative form"
379
FLExit("Velocity option spatial_discretisation/conservative_advection should be set to 0.0")
382
if (have_option(trim(velocity_option_path)//"/vertical_stabilization")) then
383
FLExit("With equation type ShallowWater you cannot use the option vertical_stabilization")
386
! boundary conditions that don't make sense:
387
do i=1, size(forbidden_bc_types)
388
if (have_option(trim(velocity_option_path)//"/boundary_conditions/type::"// &
389
trim(forbidden_bc_types(i)))) then
390
FLExit("Can't have "//trim(forbidden_bc_types(i))//" boundary condition with ShallowWater")
395
end subroutine shallow_water_equations_check_options
397
end module shallow_water_equations
added
removed
assemble/Shallow_Water_Equations.F90
