Neko 1.99.1
A portable framework for high-order spectral element flow simulations
Loading...
Searching...
No Matches
brinkman_source_term.f90
Go to the documentation of this file.
1! Copyright (c) 2023, The Neko Authors
2! All rights reserved.
3!
4! Redistribution and use in source and binary forms, with or without
5! modification, are permitted provided that the following conditions
6! are met:
7!
8! * Redistributions of source code must retain the above copyright
9! notice, this list of conditions and the following disclaimer.
10!
11! * Redistributions in binary form must reproduce the above
12! copyright notice, this list of conditions and the following
13! disclaimer in the documentation and/or other materials provided
14! with the distribution.
15!
16! * Neither the name of the authors nor the names of its
17! contributors may be used to endorse or promote products derived
18! from this software without specific prior written permission.
19!
20! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
21! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
22! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
23! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
24! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
25! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
26! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
27! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
29! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
30! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31! POSSIBILITY OF SUCH DAMAGE.
32!
35 use aabb, only : aabb_t, get_aabb
36 use coefs, only : coef_t
38 use field, only : field_t
39 use field_list, only : field_list_t
40 use math, only : cfill_mask, pwmax2
46 use file, only : file_t
47 use json_module, only : json_file, json_core, json_value
50 use tri_mesh, only : tri_mesh_t
52 use num_types, only : rp, dp
53 use point_zone, only : point_zone_t
57 use source_term, only : source_term_t
58 use utils, only : neko_error
59 use filter, only : filter_t
60 use pde_filter, only : pde_filter_t
63 use num_types, only : sp, dp
64 use time_state, only : time_state_t
65
68 use space, only: space_t, gll
69 implicit none
70 private
71
75 type, public, extends(source_term_t) :: brinkman_source_term_t
76 private
77
79 type(field_t), pointer :: indicator_unfiltered
81 type(field_t), pointer :: indicator
83 type(field_t), pointer :: brinkman
85 class(filter_t), allocatable :: filter
86 contains
88 procedure, public, pass(this) :: init => &
91 procedure, public, pass(this) :: free => brinkman_source_term_free
93 procedure, public, pass(this) :: compute_ => brinkman_source_term_compute
94
95 ! ----------------------------------------------------------------------- !
96 ! Private methods
97 procedure, pass(this) :: init_boundary_mesh
98 procedure, pass(this) :: init_point_zone
99
101
102contains
103
104 ! ========================================================================== !
105 ! Public methods
106
112 subroutine brinkman_source_term_init_from_json(this, json, fields, coef, &
113 variable_name)
114 class(brinkman_source_term_t), intent(inout) :: this
115 type(json_file), intent(inout) :: json
116 type(field_list_t), intent(in), target :: fields
117 type(coef_t), intent(in), target :: coef
118 character(len=*), intent(in) :: variable_name
119 real(kind=rp) :: start_time, end_time
120
121 character(len=:), allocatable :: filter_type
122 real(kind=rp) :: filter_radius
123 real(kind=rp), dimension(:), allocatable :: brinkman_limits
124 real(kind=rp) :: brinkman_penalty
125
126 type(json_value), pointer :: json_object_list
127 type(json_core) :: core
128
129 character(len=:), allocatable :: object_type
130 type(json_file) :: object_settings
131 integer :: n_regions
132 integer :: i
133 type(fld_file_output_t) :: output
134
135
136 ! Mandatory fields for the general source term
137 call json_get_or_default(json, "start_time", start_time, 0.0_rp)
138 call json_get_or_default(json, "end_time", end_time, huge(0.0_rp))
139
140 ! Read the options for the permeability field
141 call json_get(json, 'brinkman.limits', brinkman_limits)
142 call json_get(json, 'brinkman.penalty', brinkman_penalty)
143
144 if (size(brinkman_limits) .ne. 2) then
145 call neko_error('brinkman_limits must be a 2 element array of reals')
146 end if
147
148 call this%free()
149 call this%init_base(fields, coef, start_time, end_time)
150
151 ! ------------------------------------------------------------------------ !
152 ! Allocate the permeability and indicator field
153
154 call neko_field_registry%add_field(coef%dof, 'brinkman_indicator', .true.)
155 call neko_field_registry%add_field(coef%dof, 'brinkman_indicator_unfiltered', &
156 .true.)
157 call neko_field_registry%add_field(coef%dof, 'brinkman_permeability', &
158 .true.)
159
160 this%indicator => neko_field_registry%get_field('brinkman_indicator')
161 this%indicator_unfiltered => &
162 neko_field_registry%get_field('brinkman_indicator_unfiltered')
163 this%brinkman => neko_field_registry%get_field('brinkman_permeability')
164
165 ! ------------------------------------------------------------------------ !
166 ! Select which constructor should be called
167
168 call json%info('objects', n_children = n_regions)
169
170 do i = 1, n_regions
171 call json_extract_item(json, "objects", i, object_settings)
172 call json_get_or_default(object_settings, 'type', object_type, 'none')
173
174 select case (object_type)
175 case ('boundary_mesh')
176 call this%init_boundary_mesh(object_settings)
177 case ('point_zone')
178 call this%init_point_zone(object_settings)
179
180 case ('none')
181 call object_settings%print()
182 call neko_error('Brinkman source term objects require a region type')
183 case default
184 call neko_error('Brinkman source term unknown region type')
185 end select
186
187 end do
188
189 ! ------------------------------------------------------------------------ !
190 ! Filter the indicator field
191
192 call json_get_or_default(json, 'filter.type', filter_type, 'none')
193 select case (filter_type)
194 case ('PDE')
195 ! Initialize the unfiltered design field
196 call this%indicator_unfiltered%init(coef%dof)
197
198 ! Allocate a PDE filter
199 allocate(pde_filter_t::this%filter)
200
201 ! Initialize the filter
202 call this%filter%init(json, coef)
203
204 ! Copy the current indicator to unfiltered (essentially a rename)
205 call field_copy(this%indicator_unfiltered, this%indicator)
206
207 ! Apply the filter
208 call this%filter%apply(this%indicator, this%indicator_unfiltered)
209
210 ! Set up sampler to include the unfiltered and filtered fields
211 call output%init(sp, 'brinkman', 3)
212 call output%fields%assign_to_field(1, this%indicator_unfiltered)
213 call output%fields%assign_to_field(2, this%indicator)
214 call output%fields%assign_to_field(3, this%brinkman)
215
216 case ('none')
217 ! Set up sampler to include the unfiltered field
218 call output%init(sp, 'brinkman', 2)
219 call output%fields%assign_to_field(1, this%indicator)
220 call output%fields%assign_to_field(2, this%brinkman)
221
222 case default
223 call neko_error('Brinkman source term unknown filter type')
224 end select
225
226 ! ------------------------------------------------------------------------ !
227 ! Compute the permeability field
228
229 this%brinkman = this%indicator
230 call permeability_field(this%brinkman, &
231 brinkman_limits(1), brinkman_limits(2), brinkman_penalty)
232
233 ! Sample the Brinkman field
234 call output%sample(0.0_rp)
235
237
240 class(brinkman_source_term_t), intent(inout) :: this
241
242 nullify(this%indicator)
243 nullify(this%indicator_unfiltered)
244 nullify(this%brinkman)
245
246 if (allocated(this%filter)) then
247 call this%filter%free()
248 deallocate(this%filter)
249 end if
250 call this%free_base()
251
252 end subroutine brinkman_source_term_free
253
256 subroutine brinkman_source_term_compute(this, time)
257 class(brinkman_source_term_t), intent(inout) :: this
258 type(time_state_t), intent(in) :: time
259 type(field_t), pointer :: u, v, w, fu, fv, fw
260 integer :: n
261
262 n = this%fields%item_size(1)
263
264 u => neko_field_registry%get_field('u')
265 v => neko_field_registry%get_field('v')
266 w => neko_field_registry%get_field('w')
267
268 fu => this%fields%get(1)
269 fv => this%fields%get(2)
270 fw => this%fields%get(3)
271
272 call field_subcol3(fu, u, this%brinkman, n)
273 call field_subcol3(fv, v, this%brinkman, n)
274 call field_subcol3(fw, w, this%brinkman, n)
275
276 end subroutine brinkman_source_term_compute
277
278 ! ========================================================================== !
279 ! Private methods
280
282 subroutine init_boundary_mesh(this, json)
283 class(brinkman_source_term_t), intent(inout) :: this
284 type(json_file), intent(inout) :: json
285
286 ! Options
287 character(len=:), allocatable :: mesh_file_name
288 character(len=:), allocatable :: distance_transform
289 character(len=:), allocatable :: filter_type
290 character(len=:), allocatable :: mesh_transform
291
292 ! Read the options for the boundary mesh
293 type(file_t) :: mesh_file
294 type(tri_mesh_t) :: boundary_mesh
295 real(kind=rp) :: scalar_r
296 real(kind=dp) :: scalar_d
297 logical :: cache, cache_exist
298 character(len=:), allocatable :: cache_filename
299 type(file_t) :: cache_file
300 type(fld_file_output_t) :: cache_output
301 type(fld_file_data_t) :: cache_data
302 type(global_interpolation_t) :: global_interp
303 type(space_t) :: prev_Xh
304 type(interpolator_t) :: space_interp
305
306 ! Mesh transform options variables
307 real(kind=dp), dimension(:), allocatable :: box_min, box_max
308 logical :: keep_aspect_ratio
309 real(kind=dp), dimension(3) :: scaling
310 real(kind=dp), dimension(3) :: translation
311 type(field_t) :: temp_field
312 type(aabb_t) :: mesh_box, target_box
313 integer :: idx_p
314 character(len=LOG_SIZE) :: log_msg
315
316 ! ------------------------------------------------------------------------ !
317 ! Read the options for the boundary mesh
318
319 call json_get(json, 'name', mesh_file_name)
320 call json_get_or_default(json, 'cache', cache, .false.)
321
322 ! Settings on how to filter the design field
323 call json_get(json, 'distance_transform.type', distance_transform)
324
325 ! ------------------------------------------------------------------------ !
326 ! Check if we can load from cache
327 if (cache) then
328 call json_get(json, 'cache_file', cache_filename)
329
330 inquire(file=trim(cache_filename) // "0.nek5000", exist=cache_exist)
331 write(log_msg, '(A)') "Checking for Brinkman source term cache."
332 call neko_log%message(log_msg, neko_log_debug)
333
334 if (cache_exist) then
335 write(log_msg, '(A)') "Loading Brinkman source term from cache."
336 call neko_log%message(log_msg, neko_log_debug)
337
338 call cache_data%init()
339 call temp_field%init(this%coef%dof)
340
341 call cache_file%init(cache_filename // "0.fld")
342 call cache_file%set_counter(0)
343 call cache_file%read(cache_data)
344
345 !
346 ! Check that the data in the fld file matches the current case.
347 ! Note that this is a safeguard and there are corner cases where
348 ! two different meshes have the same dimension and same # of elements
349 ! but this should be enough to cover obvious cases.
350 !
351 if (cache_data%glb_nelv .ne. temp_field%msh%glb_nelv .or. &
352 cache_data%gdim .ne. temp_field%msh%gdim) then
353 call neko_error("The fld file must match the current mesh! " // &
354 "Use 'interpolate': 'true' to enable interpolation.")
355 end if
356
357 ! Do the space-to-space interpolation
358 call prev_xh%init(gll, cache_data%lx, cache_data%ly, cache_data%lz)
359 call space_interp%init(temp_field%Xh, prev_xh)
360 call space_interp%map_host(temp_field%x, cache_data%p%x, &
361 cache_data%nelv, temp_field%Xh)
362 call space_interp%free()
363 call prev_xh%free()
364
365 ! Synchronize to device if needed
366 if (neko_bcknd_device .eq. 1) then
367 call device_memcpy(temp_field%x, temp_field%x_d, &
368 temp_field%size(), host_to_device, sync = .true.)
369 end if
370
371 ! Update the global indicator field by max operator
372 call field_pwmax2(this%indicator, temp_field)
373
374 ! Clean up
375 call cache_data%free()
376 call temp_field%free()
377 call cache_file%free()
378 return
379 end if
380 end if
381
382 ! ------------------------------------------------------------------------ !
383 ! Load the immersed boundary mesh
384
385 call mesh_file%init(mesh_file_name)
386 call mesh_file%read(boundary_mesh)
387
388 if (boundary_mesh%nelv .eq. 0) then
389 call neko_error('No elements in the boundary mesh')
390 end if
391
392 ! ------------------------------------------------------------------------ !
393 ! Transform the mesh if specified.
394
395 call json_get_or_default(json, 'mesh_transform.type', &
396 mesh_transform, 'none')
397
398 select case (mesh_transform)
399 case ('none')
400 ! Do nothing
401 case ('bounding_box')
402 call json_get(json, 'mesh_transform.box_min', box_min)
403 call json_get(json, 'mesh_transform.box_max', box_max)
404 call json_get_or_default(json, 'mesh_transform.keep_aspect_ratio', &
405 keep_aspect_ratio, .true.)
406
407 if (size(box_min) .ne. 3 .or. size(box_max) .ne. 3) then
408 call neko_error('Case file: mesh_transform. &
409 &box_min and box_max must be 3 element arrays of reals')
410 end if
411
412 call target_box%init(box_min, box_max)
413
414 mesh_box = get_aabb(boundary_mesh)
415
416 scaling = target_box%get_diagonal() / mesh_box%get_diagonal()
417 if (keep_aspect_ratio) then
418 scaling = minval(scaling)
419 end if
420
421 translation = - scaling * mesh_box%get_min() + target_box%get_min()
422
423 do idx_p = 1, boundary_mesh%mpts
424 boundary_mesh%points(idx_p)%x = &
425 scaling * boundary_mesh%points(idx_p)%x + translation
426 end do
427
428 ! Report the transformation applied
429 write(log_msg, '(A)') "The following transformation was applied:"
430 call neko_log%message(log_msg)
431 write(log_msg, '(A, 3F12.6)') "Scaling: ", scaling
432 call neko_log%message(log_msg)
433 write(log_msg, '(A, 3F12.6)') "Translation: ", translation
434 call neko_log%message(log_msg)
435
436 case default
437 call neko_error('Unknown mesh transform')
438 end select
439
440 ! ------------------------------------------------------------------------ !
441 ! Compute the permeability field
442
443 ! Assign the signed distance field to all GLL points in the permeability
444 ! field. Initially we just run a brute force loop over all GLL points and
445 ! compute the signed distance function. This should be replaced with a
446 ! more efficient method, such as a tree search.
447
448 call temp_field%init(this%coef%dof)
449
450 ! Select how to transform the distance field to a design field
451 select case (distance_transform)
452 case ('smooth_step')
453 call json_get(json, 'distance_transform.value', scalar_d)
454 scalar_r = real(scalar_d, kind=rp)
455
456 call signed_distance_field(temp_field, boundary_mesh, scalar_d)
457 call smooth_step_field(temp_field, scalar_r, 0.0_rp)
458
459 case ('step')
460
461 call json_get(json, 'distance_transform.value', scalar_d)
462 scalar_r = real(scalar_d, kind=rp)
463
464 call signed_distance_field(temp_field, boundary_mesh, scalar_d)
465 call step_function_field(temp_field, scalar_r, 1.0_rp, 0.0_rp)
466
467 case default
468 call neko_error('Unknown distance transform')
469 end select
470
471 ! Write the field to cache
472 if (cache) then
473 write(log_msg, '(A)') "Writing Brinkman source term to cache."
474 call neko_log%message(log_msg, neko_log_debug)
475 call cache_output%init(dp, cache_filename, 1)
476 call cache_output%fields%assign_to_field(1, temp_field)
477 call cache_output%sample(0.0_rp)
478 end if
479
480 ! Update the global indicator field by max operator
481 call field_pwmax2(this%indicator, temp_field)
482
483 call temp_field%free()
484
485 end subroutine init_boundary_mesh
486
488 subroutine init_point_zone(this, json)
489 class(brinkman_source_term_t), intent(inout) :: this
490 type(json_file), intent(inout) :: json
491
492 ! Options
493 character(len=:), allocatable :: zone_name
494
495 type(field_t) :: temp_field
496 class(point_zone_t), pointer :: zone
497 integer :: i
498
499 ! ------------------------------------------------------------------------ !
500 ! Read the options for the point zone
501
502 call json_get(json, 'name', zone_name)
503
504 ! Compute the indicator field
505 call temp_field%init(this%coef%dof)
506
507 zone => neko_point_zone_registry%get_point_zone(zone_name)
508
509 if (neko_bcknd_device .eq. 1) then
510 call device_cfill_mask(temp_field%x_d, 1.0_rp, temp_field%size(), &
511 zone%mask%get_d(), zone%size)
512 else
513 call cfill_mask(temp_field%x, 1.0_rp, temp_field%size(), &
514 zone%mask%get(), zone%size)
515 end if
516
517 ! Update the global indicator field by max operator
518 call field_pwmax2(this%indicator, temp_field)
519
520 end subroutine init_point_zone
521
522end module brinkman_source_term
double real
Copy data between host and device (or device and device)
Definition device.F90:66
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Retrieves a parameter by name or throws an error.
Axis Aligned Bounding Box (aabb) implementation in Fortran.
Definition aabb.f90:71
type(aabb_t) function, public get_aabb(object, padding)
Construct the aabb of a predefined object.
Definition aabb.f90:181
Implements the brinkman_source_term_t type.
subroutine init_point_zone(this, json)
Initializes the source term from a point zone.
subroutine init_boundary_mesh(this, json)
Initializes the source term from a boundary mesh.
subroutine brinkman_source_term_init_from_json(this, json, fields, coef, variable_name)
The common constructor using a JSON object.
subroutine brinkman_source_term_free(this)
Destructor.
subroutine brinkman_source_term_compute(this, time)
Computes the source term and adds the result to fields.
Coefficients.
Definition coef.f90:34
subroutine, public device_cfill_mask(a_d, c, n, mask_d, n_mask, strm)
Fill a constant to a masked vector. .
subroutine, public device_pwmax2(a_d, b_d, n, strm)
Compute the point-wise maximum of two vectors .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:47
subroutine, public field_pwmax2(a, b, n)
Point-wise max operation .
subroutine, public field_copy(a, b, n)
Copy a vector .
subroutine, public field_subcol3(a, b, c, n)
Returns .
Defines a registry for storing solution fields.
type(field_registry_t), target, public neko_field_registry
Global field registry.
Defines a field.
Definition field.f90:34
Module for file I/O operations.
Definition file.f90:34
Filter to be applied to a scalar field.
Definition filter.f90:38
Simple module to handle fld file series. Provides an interface to the different fields sotred in a fl...
Implements fld_file_output_t.
Implements global_interpolation given a dofmap.
Routines to interpolate between different spaces.
Utilities for retrieving parameters from the case files.
Logging routines.
Definition log.f90:34
integer, parameter, public neko_log_debug
Debug log level.
Definition log.f90:78
type(log_t), public neko_log
Global log stream.
Definition log.f90:70
integer, parameter, public log_size
Definition log.f90:46
A module containing mapping functions and subroutines. These functions are used to modify fields in a...
Definition mappings.f90:36
subroutine, public permeability_field(f, k_0, k_1, q)
Apply a permeability function to a field.
Definition mappings.f90:90
subroutine, public smooth_step_field(f, edge0, edge1)
Apply a smooth step function to a field.
Definition mappings.f90:70
subroutine, public step_function_field(f, x0, value0, value1)
Apply a step function to a field.
Definition mappings.f90:109
Definition math.f90:60
subroutine, public pwmax2(a, b, n)
Point-wise maximum of two vectors .
Definition math.f90:1385
subroutine, public cfill_mask(a, c, n, mask, n_mask)
Fill a constant to a masked vector. .
Definition math.f90:397
Build configurations.
integer, parameter neko_bcknd_device
integer, parameter, public dp
Definition num_types.f90:9
integer, parameter, public sp
Definition num_types.f90:8
integer, parameter, public rp
Global precision used in computations.
Definition num_types.f90:12
Defines an output.
Definition output.f90:34
A PDE based filter.
type(point_zone_registry_t), target, public neko_point_zone_registry
Global point_zone registry.
Profiling interface.
Definition profiler.F90:34
subroutine, public profiler_start_region(name, region_id)
Started a named (name) profiler region.
Definition profiler.F90:78
subroutine, public profiler_end_region(name, region_id)
End the most recently started profiler region.
Definition profiler.F90:109
Module containing Signed Distance Functions.
subroutine, public signed_distance_field(field_data, object, max_distance)
Signed distance field.
Implements the source_term_t type and a wrapper source_term_wrapper_t.
Defines a function space.
Definition space.f90:34
integer, parameter, public gll
Definition space.f90:48
Module with things related to the simulation time.
Defines a triangular surface mesh.
Definition tri_mesh.f90:35
Utilities.
Definition utils.f90:35
Axis Aligned Bounding Box (aabb) data structure.
Definition aabb.f90:107
A Brinkman source term. The region and strength are controlled by assigning regions types and brinkma...
Coefficients defined on a given (mesh, ) tuple. Arrays use indices (i,j,k,e): element e,...
Definition coef.f90:55
field_list_t, To be able to group fields together
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Definition file.f90:55
Base abstract class for filter.
Definition filter.f90:48
A simple output saving a list of fields to a .fld file.
Implements global interpolation for arbitrary points in the domain.
Interpolation between two space::space_t.
A PDE based filter mapping $\rho \mapsto \tilde{\rho}$, see Lazarov & O. Sigmund 2010,...
Base abstract type for point zones.
Base abstract type for source terms.
The function space for the SEM solution fields.
Definition space.f90:62
A struct that contains all info about the time, expand as needed.