35
36
use sparse_matrices_fields
38
use global_parameters, only: OPTION_PATH_LEN
39
use global_parameters, only: OPTION_PATH_LEN, dt, timestep
39
40
use state_fields_module
40
41
use boundary_conditions
45
integer, dimension(:), allocatable :: sediment_boundary_condition_ids
47
public set_sediment_reentrainment, sediment_init, set_sediment_bc_id
48
public sediment_cleanup, get_sediment_bc_id
49
public get_sediment_name, get_nSediments
46
public set_sediment_reentrainment, sediment_check_options, get_n_sediment_fields, &
51
interface get_sediment_item
52
module procedure get_sediment_field, get_sediment_field_name,&
53
& get_sediment_option_string, get_sediment_option_real,&
54
& get_sediment_option_scalar_field
55
end interface get_sediment_item
55
function get_sediment_name(i)
56
!!< Return the name of sediment class from the options tree at position i-1.
57
!!< Loops over sediment should therefore be 1,nSedimentClasses in order to
58
!!< make this funciton work.
60
integer, intent(in) :: i
62
character(len=FIELD_NAME_LEN) :: get_sediment_name
64
call get_option('/material_phase[0]/sediment/class['//int2str(i-1)//']/name',get_sediment_name)
66
end function get_sediment_name
69
function get_nSediments()
70
!!< Return the number of sediment classes
72
integer :: get_nSediments
74
get_nSediments = option_count('/material_phase[0]/sediment/class')
76
end function get_nSediments
78
subroutine sediment_init()
80
integer :: ic, bc, i, stat, nSedimentClasses
81
integer :: nBoundaryConditions, nInitialConditions, existingInitCond
82
character(len=OPTION_PATH_LEN) :: option_path, temp_path
83
character(len=FIELD_NAME_LEN) :: class_name, bc_name, ic_name
85
! Set up some constants and allocate an array of sediment boundary condition
87
nSedimentClasses = option_count('/material_phase[0]/sediment/sediment_class')
88
allocate(sediment_boundary_condition_ids(nSedimentClasses))
89
sediment_boundary_condition_ids = -1
91
ewrite(2,*) "In sediment_init"
92
ewrite(2,*) "Found ",nSedimentClasses, "sediment classes"
93
! Check if we have SedimentTemplate on - if so, this is a fresh run, so
94
! construct the sediment classes from the template.
95
! If it doesn't exist, it's a checkpoint run, so no need to do this.
96
! Note that we have to also check options more carefully, as we can't
97
! carry out all possible checks in options_check (our fields haven't been
98
! constructed at that point).
99
if (have_option('/material_phase[0]/sediment/scalar_field::SedimentTemplate')) then
101
ewrite(2,*) "Initialising sediment classes from Template"
103
! For each sediment class
104
sediment_classes: do i=1,nSedimentClasses
105
! copy the template into a temp option path
106
temp_path = '/material_phase[0]/sediment/Sediment_'//int2str(i-1)//'_temp'
107
call copy_option('/material_phase[0]/sediment/scalar_field::SedimentTemplate',&
110
! now go over all the top level options and copy them into the
112
option_path='/material_phase[0]/sediment/sediment_class['//int2str(i-1)//']'
113
call get_option(trim(option_path)//'/name',class_name)
114
option_path='/material_phase[0]/sediment/sediment_class::'//trim(class_name)
115
call delete_option(trim(temp_path)//"/name")
116
call copy_option(trim(option_path)//"/name",trim(temp_path)//"/name")
118
! loop over all boundary conditions
119
nBoundaryConditions = option_count(trim(option_path)//'/boundary_conditions')
120
bc_loop: do bc=1,nBoundaryConditions
121
call get_option(trim(option_path)//'/boundary_conditions['//int2str(bc-1)//']/name',bc_name)
122
call move_option(trim(option_path)//'/boundary_conditions::'//trim(bc_name),&
123
trim(temp_path)//'/prognostic//boundary_conditions::'//trim(bc_name))
126
! loop over all initial conditions
127
nInitialConditions = option_count(trim(option_path)//'/initial_condition')
128
existingInitCond = option_count(trim(temp_path)//'/prognostic/initial_condition')
129
if (nInitialConditions > 0) then
130
do ic=1,existingInitCond
131
call delete_option(trim(temp_path)//'/prognostic/initial_condition[0]')
134
ic_loop: do ic=1,nInitialConditions
135
call get_option(trim(option_path)//'/initial_condition['//int2str(ic-1)//']/name',ic_name)
136
call copy_option(trim(option_path)//'/initial_condition::'//trim(ic_name),&
137
trim(temp_path)//'/prognostic/initial_condition::'//trim(ic_name))
140
! now do the one off fields and options
141
if (have_option(trim(option_path)//'/tensor_field::Diffusivity')) then
142
if (have_option(trim(temp_path)//'/prognostic/tensor_field::Diffusivity')) then
143
call delete_option(trim(temp_path)//'/prognostic/tensor_field::Diffusivity')
145
call move_option(trim(option_path)//'/tensor_field::Diffusivity',&
146
& trim(temp_path)//'/prognostic/tensor_field::Diffusivity')
149
! Check that we have a density somewhere
150
if (.not. have_option(trim(option_path)//'/density') .and. &
151
.not. have_option(trim(temp_path)//'/density')) then
152
FLExit("You must specify a sediment density under the Template &&
153
&& or under each sediment class")
155
! OK, so let's make sure we have the one under the sediment class
157
if (have_option(trim(option_path)//'/density')) then
158
if (have_option(trim(temp_path)//'/density')) then
159
call delete_option(trim(temp_path)//'/density')
161
call move_option(trim(option_path)//'/density',&
162
& trim(temp_path)//'/density')
165
if (have_option(trim(option_path)//'/diameter')) then
166
if (have_option(trim(temp_path)//'/diameter')) then
167
call delete_option(trim(temp_path)//'/diameter')
169
call move_option(trim(option_path)//'/diameter',&
170
& trim(temp_path)//'/diameter')
173
if (have_option(trim(option_path)//'/erodability')) then
174
if (have_option(trim(temp_path)//'/erodability')) then
175
call delete_option(trim(temp_path)//'/erodability')
177
call move_option(trim(option_path)//'/erodability',&
178
& trim(temp_path)//'/erodability')
181
if (have_option(trim(option_path)//'/critical_shear_stress')) then
182
if (have_option(trim(temp_path)//'/critical_shear_stress')) then
183
call delete_option(trim(temp_path)//'/critical_shear_stress')
185
call move_option(trim(option_path)//'/critical_shear_stress',&
186
& trim(temp_path)//'/critical_shear_stress')
189
if (have_option(trim(option_path)//'/scalar_field::SinkingVelocity')) then
190
if (have_option(trim(temp_path)//'/prognostic/scalar_field::SinkingVelocity')) then
191
call delete_option(trim(temp_path)//'/prognostic/scalar_field::SinkingVelocity')
193
call move_option(trim(option_path)//'/scalar_field::SinkingVelocity',&
194
& trim(temp_path)//'/prognostic/scalar_field::SinkingVelocity')
197
! All done, move over the temporary one to the main material phase
198
! Note that this is no longer under /sediment, so it acts like a
199
! normal scalar field from now on...
200
call move_option(temp_path,'/material_phase[0]/scalar_field::'//trim(class_name))
201
! ...however, after a restrat from checkpoint, we won't know about
202
! the sediment fields any longer, so add them to the option tree
203
call add_option('/material_phase[0]/sediment/class['//int2str(i-1)//']',stat)
204
call add_option('/material_phase[0]/sediment/class['//int2str(i-1)//']/name',stat)
205
call set_option('/material_phase[0]/sediment/class['//int2str(i-1)//']/name',trim(class_name))
207
end do sediment_classes
209
! Remove the Sediment template from the options tree
210
call delete_option('/material_phase[0]/sediment/scalar_field::SedimentTemplate')
211
! delete the old sediment_class options
212
do i=1,nSedimentClasses
213
! only ever need to delete 0 and when we remove a field, the next
214
! one will become zero
215
call delete_option('/material_phase[0]/sediment/sediment_class[0]')
219
! restart from checkpoint, so check the modified options tree
220
nSedimentClasses = option_count('/material_phase[0]/sediment/class/')
221
ewrite(2,*) "Initialising sediment classes from pre-existing run"
222
ewrite(2,*) "Found ",nSedimentClasses, "sediment classes"
223
ewrite(2,*) "They live in: "
224
do i=1,nSedimentClasses
225
call get_option('/material_phase[0]/sediment/class['//int2str(i-1)//']/name',class_name)
226
ewrite(2,*) "Sediment class ",i," ",trim(class_name)
231
end subroutine sediment_init
233
subroutine sediment_cleanup
235
deallocate(sediment_boundary_condition_ids)
237
end subroutine sediment_cleanup
239
subroutine set_sediment_bc_id(name, id)
241
character(len=FIELD_NAME_LEN), intent(in) :: name
242
integer, intent(in) :: id
244
integer :: i, nSedimentClasses
245
character(len=FIELD_NAME_LEN) :: class_name
247
nSedimentClasses = get_nSediments()
249
do i=1,nSedimentClasses
251
class_name = get_sediment_name(i)
253
if (class_name .eq. name) then
254
sediment_boundary_condition_ids(i) = id
259
end subroutine set_sediment_bc_id
261
function get_sediment_bc_id(i) result (id)
263
integer, intent(in) :: i
267
id = sediment_boundary_condition_ids(i)
269
end function get_sediment_bc_id
271
subroutine set_sediment_reentrainment(state)
273
type(state_type), intent(in):: state
275
type(scalar_field), pointer :: SedConc, erosion, bedload
276
type(vector_field), pointer :: bedShearStress
277
real :: erodibility, porosity, critical_shear_stress, shear
278
real :: erosion_flux, diameter, density, g, s
280
integer :: nNodes, i, j, nSedimentClasses
281
character(len=FIELD_NAME_LEN) :: class_name
282
character(len=OPTION_PATH_LEN) :: option_path
283
type(scalar_field) :: bedLoadSurface
284
type(vector_field) :: shearStressSurface
285
integer, dimension(:), pointer :: surface_element_list
286
type(mesh_type), pointer :: bottom_mesh
290
ewrite(1,*) "In set_sediment_bc"
292
call get_option("/timestepping/timestep", dt)
59
function get_n_sediment_fields() result (n_fields)
61
! Returns the number of sediment fields
64
n_fields = option_count('/material_phase[0]/sediment/scalar_field')
66
if (have_option('/material_phase[0]/sediment/scalar_field::SedimentBedActiveLayer&
67
&D50')) n_fields = n_fields - 1
68
if (have_option('/material_phase[0]/sediment/scalar_field::SedimentBedActiveLayer&
69
&Sigma')) n_fields = n_fields - 1
70
if (have_option('/material_phase[0]/sediment/scalar_field::ZeroSedimentConcentrat&
71
&ionViscosity')) n_fields = n_fields - 1
73
end function get_n_sediment_fields
75
subroutine get_sediment_field(state, i_field, item, stat)
77
! Returns sediment field string option
78
type(state_type), intent(in) :: state
79
integer, intent(in) :: i_field
80
type(scalar_field), pointer, intent(out) :: item
81
integer, intent(out), optional :: stat
82
character(len=FIELD_NAME_LEN) :: name
84
call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
86
item => extract_scalar_field(state, trim(name), stat)
88
end subroutine get_sediment_field
90
subroutine get_sediment_field_name(state, i_field, item, stat)
92
! Returns sediment field string option
93
type(state_type), intent(in) :: state
94
integer, intent(in) :: i_field
95
character(len=FIELD_NAME_LEN), intent(out) :: item
96
integer, intent(out), optional :: stat
98
call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
99
& 1)//']/name', item, stat=stat)
101
end subroutine get_sediment_field_name
103
subroutine get_sediment_option_string(state, i_field, option, item, stat)
105
! Returns sediment field string option
106
type(state_type), intent(in) :: state
107
integer, intent(in) :: i_field
108
character(len=*), intent(in) :: option
109
character(len=FIELD_NAME_LEN), intent(out) :: item
110
integer, intent(out), optional :: stat
112
call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
113
& 1)//']/prognostic/'//option, item, stat=stat)
115
end subroutine get_sediment_option_string
117
subroutine get_sediment_option_real(state, i_field, option, item, stat, default)
119
! Returns sediment field real option
120
type(state_type), intent(in) :: state
121
integer, intent(in) :: i_field
122
character(len=*), intent(in) :: option
123
real, intent(out) :: item
124
integer, intent(out), optional :: stat
125
real, intent(in), optional :: default
127
call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
128
& 1)//']/prognostic/'//option, item, stat = stat, default = default)
130
end subroutine get_sediment_option_real
132
subroutine get_sediment_option_scalar_field(state, i_field, option, item, stat)
134
! Returns sediment field related scalar field
135
type(state_type), intent(in) :: state
136
integer, intent(in) :: i_field
137
type(scalar_field), pointer, intent(out) :: item
138
character(len=*), intent(in) :: option
139
integer, intent(out), optional :: stat
141
character(len=FIELD_NAME_LEN) :: field_name
143
call get_option(trim(state%option_path)//'/sediment/scalar_field['//int2str(i_field -&
144
& 1)//']/name', field_name)
145
item => extract_scalar_field(state, trim(field_name)//option, stat)
147
end subroutine get_sediment_option_scalar_field
149
subroutine set_sediment_reentrainment(state)
151
type(state_type), intent(in) :: state
152
type(scalar_field), pointer :: sediment_field
153
integer :: i_field, i_bc, n_bc
154
character(len = FIELD_NAME_LEN) :: bc_name, bc_type
155
character(len = OPTION_PATH_LEN) :: bc_path, bc_path_
157
ewrite(1,*) "In set_sediment_reentrainment"
159
sediment_fields: do i_field = 1, get_n_sediment_fields()
161
! extract sediment field from state
162
call get_sediment_item(state, i_field, sediment_field)
164
! get boundary condition path and number of boundary conditions
165
bc_path = trim(sediment_field%option_path)//'/prognostic/boundary_conditions'
166
n_bc = option_count(trim(bc_path))
168
! Loop over boundary conditions for field
169
boundary_conditions: do i_bc = 0, n_bc - 1
171
! Get name and type of boundary condition
172
bc_path_ = trim(bc_path)//"["//int2str(i_bc)//"]"
173
call get_option(trim(bc_path_)//"/name", bc_name)
174
call get_option(trim(bc_path_)//"/type[0]/name", bc_type)
176
! skip if this is not a reentrainment boundary
177
if (.not. (trim(bc_type) .eq. "sediment_reentrainment")) then
178
cycle boundary_conditions
181
ewrite(1,*) "Setting reentrainment boundary condition "//trim(bc_name)//" for fi&
182
&eld: "//trim(sediment_field%name)
184
call set_reentrainment_bc(state, sediment_field, bc_name, bc_path_, i_field)
186
end do boundary_conditions
188
end do sediment_fields
190
end subroutine set_sediment_reentrainment
192
subroutine set_reentrainment_bc(state, sediment_field, bc_name, bc_path, i_field)
194
type(state_type), intent(in) :: state
195
type(scalar_field), intent(in), pointer :: sediment_field
196
type(scalar_field), pointer :: reentrainment, bedload, sink_U, d50,&
197
& sigma, volume_fraction
198
type(scalar_field) :: masslump, bedload_remap
199
type(tensor_field), pointer :: viscosity_pointer
200
type(tensor_field), target :: viscosity
201
type(vector_field), pointer :: x, shear_stress
202
type(mesh_type), pointer :: surface_mesh
203
character(len = FIELD_NAME_LEN) :: bc_name, bc_path, algorithm
204
integer :: stat, i_ele, i_field, i_node, i, j
205
integer, dimension(:), pointer :: surface_element_list
206
real, dimension(2,2) :: algorithm_viscosity
208
! get boundary condition field and zero
209
reentrainment => extract_surface_field(sediment_field, bc_name, 'value')
210
call set(reentrainment, 0.0)
212
! get boundary condition info
213
call get_boundary_condition(sediment_field, name=bc_name,&
214
& surface_mesh=surface_mesh, surface_element_list=surface_element_list)
217
call get_sediment_item(state, i_field, 'Bedload', bedload)
219
! get volume fraction
220
call get_sediment_item(state, i_field, 'BedloadVolumeFraction', volume_fraction)
222
! get sinking velocity
223
call get_sediment_item(state, i_field, 'SinkingVelocity', sink_U)
226
d50 => extract_scalar_field(state, 'SedimentBedActiveLayerD50', stat)
229
sigma => extract_scalar_field(state, 'SedimentBedActiveLayerSigma', stat)
231
call allocate(viscosity, sediment_field%mesh, "Viscosity")
233
call get_option(trim(bc_path)//"/type[0]/viscosity", algorithm_viscosity(1,1), stat&
238
algorithm_viscosity(i,j) = algorithm_viscosity(1,1)
242
call set(viscosity, algorithm_viscosity)
243
viscosity_pointer => viscosity
245
viscosity_pointer => extract_tensor_field(state, "Viscosity", stat)
247
FLExit("A viscosity must be specified to calculate reentrainment")
252
shear_stress => extract_vector_field(state, "BedShearStress", stat)
254
FLExit("A bed shear stress must be specified to calculate reentrainment")
257
! get coordinate field
258
x => extract_vector_field(state, "Coordinate")
260
if (continuity(surface_mesh)>=0) then
261
! For continuous fields we need a global lumped mass. For dg we'll
262
! do the mass inversion on a per face basis inside the element loop.
263
! Continuity must be the same for all bedload meshes
264
call allocate(masslump, surface_mesh, "SurfaceMassLump")
268
call get_option(trim(bc_path)//"/type[0]/algorithm", algorithm)
269
! loop through elements in surface field and calculate reentrainment
270
elements: do i_ele = 1, element_count(reentrainment)
272
select case(trim(algorithm))
274
call assemble_generic_reentrainment_ele(state, i_field, i_ele, reentrainment,&
275
& shear_stress, surface_element_list, x, masslump, sink_U, volume_fraction)
277
call assemble_garcia_1991_reentrainment_ele(state, i_field, i_ele, reentrainment,&
278
& x, masslump, surface_element_list, viscosity_pointer,&
279
& shear_stress, d50, sink_U, sigma, volume_fraction)
281
FLExit("A valid re-entrainment algorithm must be selected")
286
! invert global lumped mass for continuous fields
287
if(continuity(surface_mesh)>=0) then
288
where (masslump%val/=0.0)
289
masslump%val=1./masslump%val
291
call scale(reentrainment, masslump)
292
call deallocate(masslump)
295
! check bound of entrainment so that it does not exceed the available sediment in the
296
! bed and is larger than zero.
297
call allocate(bedload_remap, surface_mesh, name="bedload_remap")
298
call remap_field_to_surface(bedload, bedload_remap, surface_element_list)
299
nodes: do i_node = 1, node_count(reentrainment)
301
call set(reentrainment, i_node, min(max(node_val(reentrainment, i_node), 0.0),&
302
& node_val(bedload_remap, i_node)/dt))
304
call set(reentrainment, i_node, max(node_val(reentrainment, i_node), 0.0))
308
ewrite_minmax(bedload)
309
ewrite_minmax(reentrainment)
311
call deallocate(bedload_remap)
312
call deallocate(viscosity)
314
end subroutine set_reentrainment_bc
316
subroutine assemble_garcia_1991_reentrainment_ele(state, i_field, i_ele, reentrainment,&
317
& x, masslump, surface_element_list, viscosity, shear_stress, d50,&
318
& sink_U, sigma, volume_fraction)
320
type(state_type), intent(in) :: state
321
integer, intent(in) :: i_ele, i_field
322
type(tensor_field), pointer, intent(in) :: viscosity
323
type(vector_field), intent(in) :: x, shear_stress
324
type(scalar_field), intent(inout) :: masslump
325
type(scalar_field), pointer, intent(inout) :: reentrainment
326
type(scalar_field), pointer, intent(in) :: d50, sink_U, sigma,&
328
integer, dimension(:), pointer, intent(in) :: surface_element_list
329
type(element_type), pointer :: shape
330
integer, dimension(:), pointer :: ele
331
real, dimension(ele_ngi(reentrainment, i_ele)) :: detwei
332
real, dimension(ele_loc(reentrainment, i_ele), &
333
& ele_loc(reentrainment, i_ele)) :: invmass
334
real :: A, R, d, g, density
335
real, dimension(ele_ngi(reentrainment, i_ele)) :: R_p, u_star, Z
336
real, dimension(ele_loc(reentrainment, i_ele)) :: E
337
real, dimension(ele_ngi(reentrainment, i_ele)) :: shear, lambda_m
338
real, dimension(shear_stress%dim, &
339
& ele_ngi(reentrainment, i_ele)) :: shear_quad
344
ele => ele_nodes(reentrainment, i_ele)
345
shape => ele_shape(reentrainment, i_ele)
347
call transform_facet_to_physical(x, surface_element_list(i_ele), detwei)
349
if(continuity(reentrainment)>=0) then
350
call addto(masslump, ele, &
351
sum(shape_shape(shape, shape, detwei), 1))
353
! In the DG case we will apply the inverse mass locally.
354
invmass=inverse(shape_shape(shape, shape, detwei))
357
! calculate particle Reynolds number
358
call get_sediment_item(state, i_field, 'submerged_specific_gravity', R)
359
call get_sediment_item(state, i_field, 'diameter', d)
295
360
call get_option("/physical_parameters/gravity/magnitude", g)
296
viscosity = 1. / 1000.
300
nSedimentClasses = get_nSediments()
302
do i=1,nSedimentClasses
304
if (sediment_boundary_condition_ids(i) .eq. -1) then
308
class_name = get_sediment_name(i)
309
SedConc => extract_scalar_field(state, trim(class_name))
310
option_path = SedConc%option_path
311
bedload => extract_scalar_field(state,"SedimentFlux"//trim(class_name))
313
call get_option(trim(option_path)//"/porosity", porosity, default=0.3)
314
call get_option(trim(option_path)//"/name", class_name)
316
if (.not. alloced) then
317
bedShearStress => extract_vector_field(state, "BedShearStress")
318
call get_boundary_condition(SedConc, sediment_boundary_condition_ids(i), &
319
surface_mesh=bottom_mesh,surface_element_list=surface_element_list)
320
!call create_surface_mesh(bottom_mesh, surface_nodes, bedShearStress%mesh, surface_element_list, 'ErosionBed')
321
call allocate(bedLoadSurface,bottom_mesh, name="bedLoadSurface")
322
call allocate(shearStressSurface,bedShearStress%dim,bottom_mesh, name="shearStressSurface")
323
call remap_field_to_surface(bedShearStress, shearStressSurface, &
324
surface_element_list)
326
nNodes = node_count(bedLoadSurface)
329
call get_option(trim(option_path)//"/erodability", erodibility, default=1.0)
330
if (have_option(trim(option_path)//"/critical_shear_stress")) then
331
call get_option(trim(option_path)//"/critical_shear_stress", critical_shear_stress)
333
if (have_option(trim(option_path)//"/diameter")) then
334
call get_option(trim(option_path)//"/diameter",diameter)
336
FLExit("You need to either specify a critical shear stress or a &&
337
&& sediment diameter")
339
call get_option(trim(option_path)//"/density",density)
340
! calc critical shear stress
342
!S_star = sqrt((s-1)*g*diameter**3)/viscosity
343
!critical_shear_stress = 0.105*S_star**(-0.13) + &
344
! 0.045*exp(-35*S_star**(-0.59))
345
! estimate of critical shear stress assuming grains larger than
346
! 10 microns and constant viscosity - note the conversion to mm!
347
critical_shear_stress = 0.041 * (s-1) * 1024. * g * (diameter/1000.)
350
call remap_field_to_surface(bedload, bedLoadSurface, &
351
surface_element_list)
352
erosion => extract_surface_field(SedConc,sediment_boundary_condition_ids(i),"value")
354
! we only need to add to the source the erosion of sediment from the
355
! bedload into the neumann BC term
357
! The source depends on the erodability of that sediment
358
! (usually 1 unless you want to do something like mud which sticks),
359
! the porosity in the bedload and the bed shear stress.
361
! Each sediment class has a critical shear stress, which if exceeded
362
! by the bed shear stress, sediment is placed into suspension
364
! loop over nodes in bottom surface
367
shear = norm2(node_val(ShearStressSurface, j))
368
! critical stress is either given by user (optional) or calculated
369
! using Shield's formula (depends on grain size and density and
370
! (vertical) viscosity)
371
erosion_flux = erodibility*(1-porosity)*((shear - critical_shear_stress) / critical_shear_stress)
373
if (erosion_flux < 0) then
376
! A limit is placed depending on how much of that sediment is in the
378
if (erosion_flux*dt > node_val(bedLoadSurface,j)) then
379
erosion_flux = node_val(bedLoadSurface,j)/dt
382
call set(erosion,j,erosion_flux)
388
! leave it up to the standard solve to add this term in as a neumann BC.
390
call deallocate(bedLoadSurface)
391
call deallocate(shearStressSurface)
393
end subroutine set_sediment_reentrainment
361
! VISCOSITY ASSUMED TO BE ISOTROPIC - maybe should be in normal direction to surface
362
where (face_val_at_quad(viscosity, surface_element_list(i_ele), 1, 1) > 0.0)
363
R_p = sqrt(R*g*d**3.0)/face_val_at_quad(viscosity, surface_element_list(i_ele), 1, 1)
368
! calculate u_star (shear velocity)
369
call get_option(trim(shear_stress%option_path)//"/diagnostic/density", density)
370
! calculate magnitude of shear stress at quadrature points
371
shear_quad = face_val_at_quad(shear_stress, surface_element_list(i_ele))
372
do i_gi = 1, ele_ngi(reentrainment, i_ele)
373
shear(i_gi) = norm2(shear_quad(:, i_gi))
375
u_star = sqrt(shear/density)
378
lambda_m = 1.0 - 0.288 * face_val_at_quad(sigma, surface_element_list(i_ele))
381
where (face_val_at_quad(d50, surface_element_list(i_ele)) > 0.0)
382
Z = lambda_m * u_star/face_val_at_quad(sink_U, surface_element_list(i_ele)) * R_p&
383
&**0.6 * (d / face_val_at_quad(d50, surface_element_list(i_ele)))**0.2
388
! calculate reentrainment F*v_s*E
389
E = shape_rhs(shape, face_val_at_quad(volume_fraction, surface_element_list(i_ele)) *&
390
& face_val_at_quad(sink_U, surface_element_list(i_ele)) * A*Z**5 / (1 + A*Z**5&
393
if(continuity(reentrainment)<0) then
395
E = matmul(invmass, E)
398
call addto(reentrainment, ele, E)
400
end subroutine assemble_garcia_1991_reentrainment_ele
402
subroutine assemble_generic_reentrainment_ele(state, i_field, i_ele, reentrainment,&
403
& shear_stress, surface_element_list, x, masslump, sink_U, volume_fraction)
405
type(state_type), intent(in) :: state
406
integer, intent(in) :: i_ele, i_field
407
type(scalar_field), pointer, intent(inout) :: reentrainment
408
type(vector_field), intent(in) :: x, shear_stress
409
integer, dimension(:), pointer, intent(in) :: surface_element_list
410
type(scalar_field), intent(inout) :: masslump
411
type(scalar_field), pointer, intent(in) :: sink_U, volume_fraction
412
type(element_type), pointer :: shape
413
integer, dimension(:), pointer :: ele
414
real, dimension(ele_ngi(reentrainment, i_ele)) :: detwei
415
real, dimension(ele_loc(reentrainment, i_ele), &
416
& ele_loc(reentrainment, i_ele)) :: invmass
417
integer, dimension(2) :: stat
418
real :: shear_crit, d, R, g, erod,&
419
& poro, density, rho_0
420
real, dimension(ele_loc(reentrainment, i_ele)) :: E
421
real, dimension(ele_ngi(reentrainment, i_ele)) :: shear
422
real, dimension(shear_stress%dim, &
423
& ele_ngi(reentrainment, i_ele)) :: shear_quad
427
call get_sediment_item(state, i_field, 'critical_shear_stress', shear_crit, stat(1))
428
call get_sediment_item(state, i_field, 'diameter', d, stat(2))
429
! non-dimensionalise shear stress
430
call get_option('/material_phase::'//trim(state%name)//'/equation_of_state/fluids/line&
431
&ar/reference_density', rho_0)
432
! get or calculate critical shear stress
433
if (.not.any(stat .eq. 0)) then
434
FLExit("You need to either specify a critical shear stress or a &
435
&sediment diameter to use the generic formula for reentrainment")
436
else if (stat(1) /= 0) then
437
! estimate of critical shear stress assuming grains larger than
438
! 10 microns and constant viscosity
439
! critical stress is either given by user (optional) or calculated
440
! using Shield's formula (depends on grain size and density and
441
! (vertical) viscosity)
442
call get_sediment_item(state, i_field, 'submerged_specific_gravity', R)
443
call get_option("/physical_parameters/gravity/magnitude", g)
444
shear_crit = 0.041 * R * rho_0 * g * d
447
! calculate eroded sediment flux and set reentrainment BC
448
! we only need to add to the source the reentrainment of sediment from the
449
! bedload into the neumann BC term
451
! The source depends on the erodability of that sediment
452
! (usually 1 unless you want to do something like mud which sticks),
453
! the porosity in the bedload and the bed shear stress.
455
! Each sediment class has a critical shear stress, which if exceeded
456
! by the bed shear stress, sediment is placed into suspension
458
! loop over nodes in bottom surface
459
call get_sediment_item(state, i_field, 'erodability', erod, default=1.0)
460
call get_sediment_item(state, i_field, 'bed_porosity', poro, default=0.3)
462
ele => ele_nodes(reentrainment, i_ele)
463
shape => ele_shape(reentrainment, i_ele)
465
call transform_facet_to_physical(x, surface_element_list(i_ele), detwei)
467
if(continuity(reentrainment)>=0) then
468
call addto(masslump, ele, &
469
sum(shape_shape(shape, shape, detwei), 1))
471
! In the DG case we will apply the inverse mass locally.
472
invmass=inverse(shape_shape(shape, shape, detwei))
475
! calculate magnitude of shear stress at quadrature points
476
shear_quad = face_val_at_quad(shear_stress, surface_element_list(i_ele))
477
do i_gi = 1, ele_ngi(reentrainment, i_ele)
478
shear(i_gi) = norm2(shear_quad(:, i_gi))
480
! non-dimensionalise shear stress
481
call get_option(trim(shear_stress%option_path)//"/diagnostic/density", density)
482
shear = shear / density
484
! calculate reentrainment F*vs*E
485
E = shape_rhs(shape, face_val_at_quad(volume_fraction, surface_element_list(i_ele)) *&
486
& face_val_at_quad(sink_U, surface_element_list(i_ele)) * erod * (1-poro) *&
487
& (shear - shear_crit)/shear_crit * detwei)
489
if(continuity(reentrainment)<0) then
491
E = matmul(invmass, E)
494
call addto(reentrainment, ele, E)
496
end subroutine assemble_generic_reentrainment_ele
498
subroutine sediment_check_options
500
character(len=FIELD_NAME_LEN) :: field_mesh, sediment_mesh, bc_type
501
character(len=OPTION_PATH_LEN) :: field_option_path
502
character(len=10) :: type
503
integer :: i_field, i_bc, i_bc_surf, i_bedload_surf,&
504
& n_sediment_fields, nbcs
505
integer, dimension(2) :: bc_surface_id_count, bedload_surface_id_count
506
integer, dimension(:), allocatable :: bc_surface_ids, bedload_surface_ids
508
if (have_option('/material_phase[0]/sediment/')) then
510
ewrite(1,*) 'Checking sediment model options'
512
n_sediment_fields = get_n_sediment_fields()
514
call get_option('/material_phase[0]/sediment/scalar_field[0]/prognostic/mesh[0]/name', sediment_mesh)
516
sediment_fields: do i_field=1,n_sediment_fields
518
field_option_path = '/material_phase[0]/sediment/scalar_field['//int2str(i_field - 1)//']/pro&
521
! check sinking velocity is specified for every sediment field
522
if (.not.(have_option(trim(field_option_path)//'/scalar_field::SinkingVelocity')))&
524
FLExit("You must specify a sinking velocity for every sediment field")
527
! check all sediment fields are on the same mesh
528
call get_option(trim(field_option_path)//'/mesh[0]/name', field_mesh)
529
if (.not.(trim(field_mesh) .eq. trim(sediment_mesh))) then
530
FLExit("All sediment fields must be on the same mesh")
533
! check re-entrainment options
534
! get boundary condition path and number of boundary conditions
535
nbcs=option_count(trim(field_option_path)//'/boundary_conditions')
536
! Loop over boundary conditions for field
537
boundary_conditions: do i_bc=0, nbcs-1
539
! get name and type of boundary condition
540
call get_option(trim(field_option_path)//'/boundary_conditions['//int2str(i_bc)//&
541
&']/type[0]/name', bc_type)
543
! check whether this is a reentrainment boundary
544
if (.not. (trim(bc_type) .eq. "sediment_reentrainment")) then
545
cycle boundary_conditions
548
! check a 'BedShearStress' field exists
549
if (.not.(have_option('/material_phase[0]/vector_field::BedShearStress'))) then
550
FLExit("Reentrainment boundary condition requires a BedShearStress field")
553
! warn if bedload field is prescribed
555
if (have_option(trim(field_option_path)//'/scalar_field::Bedload/prescribed')) then
556
ewrite(0,*) 'WARNING: Bedload field is prescribed'
560
! check boundary id's are the same for re-entrainment and bedload
562
! get bedload surface ids
563
bedload_surface_id_count=option_shape(trim(field_option_path)// &
564
'/scalar_field::Bedload/'//type//'/surface_ids')
565
allocate(bedload_surface_ids(bedload_surface_id_count(1)))
566
call get_option(trim(field_option_path)// &
567
'/scalar_field::Bedload/'//type//'/surface_ids', bedload_surface_ids)
569
! get reentrainment surface ids
570
bc_surface_id_count=option_shape(trim(field_option_path)//'/boundary_conditions['&
571
&//int2str(i_bc)//']/surface_ids')
572
allocate(bc_surface_ids(bc_surface_id_count(1)))
573
call get_option(trim(field_option_path)// &
574
'/boundary_conditions['//int2str(i_bc)//']/surface_ids', bc_surface_ids)
576
bc_surface_id: do i_bc_surf=1, bc_surface_id_count(1)
578
bedload_surface_id: do i_bedload_surf=1, bedload_surface_id_count(1)
580
if (bc_surface_ids(i_bc_surf) .eq. bedload_surface_ids(i_bedload_surf)) then
584
end do bedload_surface_id
586
FLExit("Reentrainment boundary condition is specified on a surface with no bedload")
590
deallocate(bc_surface_ids)
591
deallocate(bedload_surface_ids)
593
end do boundary_conditions
595
end do sediment_fields
597
ewrite(1,*) 'Sediment model options check complete'
601
end subroutine sediment_check_options
395
603
end module sediment
added
removed
sediments/Sediment.F90
