Loading [MathJax]/extensions/tex2jax.js
Neko 0.9.99
A portable framework for high-order spectral element flow simulations
All Classes Namespaces Files Functions Variables Typedefs Enumerator Macros Pages
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
39 use field, only: field_t
40 use field_list, only: field_list_t
44 use file, only: file_t
45 use json_module, only: json_file, json_core, json_value
47 use logger, only: neko_log, log_size
48 use math, only: pwmax, cfill_mask
49 use tri_mesh, only: tri_mesh_t
51 use num_types, only: rp, dp
52 use point_zone, only: point_zone_t
56 use source_term, only: source_term_t
57 use utils, only: neko_error
58 use filter, only: filter_t
59 use pde_filter, only: pde_filter_t
61 use num_types, only: sp
62 implicit none
63 private
64
68 type, public, extends(source_term_t) :: brinkman_source_term_t
69 private
70
72 type(field_t) :: indicator_unfiltered
74 type(field_t) :: indicator
76 type(field_t) :: brinkman
78 class(filter_t), allocatable :: filter
79 contains
81 procedure, public, pass(this) :: init => &
84 procedure, public, pass(this) :: free => brinkman_source_term_free
86 procedure, public, pass(this) :: compute_ => brinkman_source_term_compute
87
88 ! ----------------------------------------------------------------------- !
89 ! Private methods
90 procedure, pass(this) :: init_boundary_mesh
91 procedure, pass(this) :: init_point_zone
92
94
95contains
96
97 ! ========================================================================== !
98 ! Public methods
99
104 subroutine brinkman_source_term_init_from_json(this, json, fields, coef)
105 class(brinkman_source_term_t), intent(inout) :: this
106 type(json_file), intent(inout) :: json
107 type(field_list_t), intent(in), target :: fields
108 type(coef_t), intent(in), target :: coef
109 real(kind=rp) :: start_time, end_time
110
111 character(len=:), allocatable :: filter_type
112 real(kind=rp) :: filter_radius
113 real(kind=rp), dimension(:), allocatable :: brinkman_limits
114 real(kind=rp) :: brinkman_penalty
115
116 type(json_value), pointer :: json_object_list
117 type(json_core) :: core
118
119 character(len=:), allocatable :: object_type
120 type(json_file) :: object_settings
121 integer :: n_regions
122 integer :: i
123 type(fld_file_output_t) :: output
124
125
126 ! Mandatory fields for the general source term
127 call json_get_or_default(json, "start_time", start_time, 0.0_rp)
128 call json_get_or_default(json, "end_time", end_time, huge(0.0_rp))
129
130 ! Read the options for the permeability field
131 call json_get(json, 'brinkman.limits', brinkman_limits)
132 call json_get(json, 'brinkman.penalty', brinkman_penalty)
133
134 if (size(brinkman_limits) .ne. 2) then
135 call neko_error('brinkman_limits must be a 2 element array of reals')
136 end if
137
138 call this%free()
139 call this%init_base(fields, coef, start_time, end_time)
140
141 ! ------------------------------------------------------------------------ !
142 ! Allocate the permeability and indicator field
143
144 if (neko_field_registry%field_exists('brinkman_indicator') &
145 .or. neko_field_registry%field_exists('brinkman')) then
146 call neko_error('Brinkman field already exists.')
147 end if
148
149 call this%indicator%init(coef%dof)
150 call this%brinkman%init(coef%dof)
151
152 ! ------------------------------------------------------------------------ !
153 ! Select which constructor should be called
154
155 call json%get('objects', json_object_list)
156 call json%info('objects', n_children = n_regions)
157 call json%get_core(core)
158
159 do i = 1, n_regions
160 call json_extract_item(core, json_object_list, i, object_settings)
161 call json_get_or_default(object_settings, 'type', object_type, 'none')
162
163 select case (object_type)
164 case ('boundary_mesh')
165 call this%init_boundary_mesh(object_settings)
166 case ('point_zone')
167 call this%init_point_zone(object_settings)
168
169 case ('none')
170 call object_settings%print()
171 call neko_error('Brinkman source term objects require a region type')
172 case default
173 call neko_error('Brinkman source term unknown region type')
174 end select
175
176 end do
177
178 ! ------------------------------------------------------------------------ !
179 ! Filter the indicator field
180
181 call json_get_or_default(json, 'filter.type', filter_type, 'none')
182 select case (filter_type)
183 case ('PDE')
184 ! Initialize the unfiltered design field
185 call this%indicator_unfiltered%init(coef%dof)
186
187 ! Allocate a PDE filter
188 allocate(pde_filter_t::this%filter)
189
190 ! Initialize the filter
191 call this%filter%init(json, coef)
192
193 ! Copy the current indicator to unfiltered (essentially a rename)
194 call field_copy(this%indicator_unfiltered, this%indicator)
195
196 ! Apply the filter
197 call this%filter%apply(this%indicator, this%indicator_unfiltered)
198
199 ! Set up sampler to include the unfiltered and filtered fields
200 call output%init(sp, 'brinkman', 3)
201 call output%fields%assign_to_field(1, this%indicator_unfiltered)
202 call output%fields%assign_to_field(2, this%indicator)
203 call output%fields%assign_to_field(3, this%brinkman)
204
205 case ('none')
206 ! Set up sampler to include the unfiltered field
207 call output%init(sp, 'brinkman', 2)
208 call output%fields%assign_to_field(1, this%indicator)
209 call output%fields%assign_to_field(2, this%brinkman)
210
211 case default
212 call neko_error('Brinkman source term unknown filter type')
213 end select
214
215 ! ------------------------------------------------------------------------ !
216 ! Compute the permeability field
217
218 this%brinkman = this%indicator
219 call permeability_field(this%brinkman, &
220 brinkman_limits(1), brinkman_limits(2), brinkman_penalty)
221
222 ! Sample the Brinkman field
223 call output%sample(0.0_rp)
224
226
229 class(brinkman_source_term_t), intent(inout) :: this
230
231 call this%indicator%free()
232 call this%brinkman%free()
233 call this%free_base()
234
235 end subroutine brinkman_source_term_free
236
240 subroutine brinkman_source_term_compute(this, t, tstep)
241 class(brinkman_source_term_t), intent(inout) :: this
242 real(kind=rp), intent(in) :: t
243 integer, intent(in) :: tstep
244 type(field_t), pointer :: u, v, w, fu, fv, fw
245 integer :: n
246
247 n = this%fields%item_size(1)
248
249 u => neko_field_registry%get_field('u')
250 v => neko_field_registry%get_field('v')
251 w => neko_field_registry%get_field('w')
252
253 fu => this%fields%get(1)
254 fv => this%fields%get(2)
255 fw => this%fields%get(3)
256
257 call field_subcol3(fu, u, this%brinkman, n)
258 call field_subcol3(fv, v, this%brinkman, n)
259 call field_subcol3(fw, w, this%brinkman, n)
260
261 end subroutine brinkman_source_term_compute
262
263 ! ========================================================================== !
264 ! Private methods
265
267 subroutine init_boundary_mesh(this, json)
268 class(brinkman_source_term_t), intent(inout) :: this
269 type(json_file), intent(inout) :: json
270
271 ! Options
272 character(len=:), allocatable :: mesh_file_name
273 character(len=:), allocatable :: distance_transform
274 character(len=:), allocatable :: filter_type
275 character(len=:), allocatable :: mesh_transform
276
277 ! Read the options for the boundary mesh
278 type(file_t) :: mesh_file
279 type(tri_mesh_t) :: boundary_mesh
280 real(kind=rp) :: scalar_r
281 real(kind=dp) :: scalar_d
282
283 ! Mesh transform options variables
284 real(kind=dp), dimension(:), allocatable :: box_min, box_max
285 logical :: keep_aspect_ratio
286 real(kind=dp), dimension(3) :: scaling
287 real(kind=dp), dimension(3) :: translation
288 type(field_t) :: temp_field
289 type(aabb_t) :: mesh_box, target_box
290 integer :: idx_p
291 character(len=LOG_SIZE) :: log_msg
292
293 ! ------------------------------------------------------------------------ !
294 ! Read the options for the boundary mesh
295
296 call json_get(json, 'name', mesh_file_name)
297
298 ! Settings on how to filter the design field
299 call json_get(json, 'distance_transform.type', distance_transform)
300
301 ! ------------------------------------------------------------------------ !
302 ! Load the immersed boundary mesh
303
304 mesh_file = file_t(mesh_file_name)
305 call mesh_file%read(boundary_mesh)
306
307 if (boundary_mesh%nelv .eq. 0) then
308 call neko_error('No elements in the boundary mesh')
309 end if
310
311 ! ------------------------------------------------------------------------ !
312 ! Transform the mesh if specified.
313
314 call json_get_or_default(json, 'mesh_transform.type', &
315 mesh_transform, 'none')
316
317 select case (mesh_transform)
318 case ('none')
319 ! Do nothing
320 case ('bounding_box')
321 call json_get(json, 'mesh_transform.box_min', box_min)
322 call json_get(json, 'mesh_transform.box_max', box_max)
323 call json_get_or_default(json, 'mesh_transform.keep_aspect_ratio', &
324 keep_aspect_ratio, .true.)
325
326 if (size(box_min) .ne. 3 .or. size(box_max) .ne. 3) then
327 call neko_error('Case file: mesh_transform. &
328 &box_min and box_max must be 3 element arrays of reals')
329 end if
330
331 call target_box%init(box_min, box_max)
332
333 mesh_box = get_aabb(boundary_mesh)
334
335 scaling = target_box%get_diagonal() / mesh_box%get_diagonal()
336 if (keep_aspect_ratio) then
337 scaling = minval(scaling)
338 end if
339
340 translation = - scaling * mesh_box%get_min() + target_box%get_min()
341
342 do idx_p = 1, boundary_mesh%mpts
343 boundary_mesh%points(idx_p)%x = &
344 scaling * boundary_mesh%points(idx_p)%x + translation
345 end do
346
347 ! Report the transformation applied
348 write(log_msg, '(A)') "The following transformation was applied:"
349 call neko_log%message(log_msg)
350 write(log_msg, '(A, 3F12.6)') "Scaling: ", scaling
351 call neko_log%message(log_msg)
352 write(log_msg, '(A, 3F12.6)') "Translation: ", translation
353 call neko_log%message(log_msg)
354
355 case default
356 call neko_error('Unknown mesh transform')
357 end select
358
359 ! ------------------------------------------------------------------------ !
360 ! Compute the permeability field
361
362 ! Assign the signed distance field to all GLL points in the permeability
363 ! field. Initally we just run a brute force loop over all GLL points and
364 ! compute the signed distance function. This should be replaced with a
365 ! more efficient method, such as a tree search.
366
367 call temp_field%init(this%coef%dof)
368
369 ! Select how to transform the distance field to a design field
370 select case (distance_transform)
371 case ('smooth_step')
372 call json_get(json, 'distance_transform.value', scalar_d)
373 scalar_r = real(scalar_d, kind=rp)
374
375 call signed_distance_field(temp_field, boundary_mesh, scalar_d)
376 call smooth_step_field(temp_field, scalar_r, 0.0_rp)
377
378 case ('step')
379
380 call json_get(json, 'distance_transform.value', scalar_d)
381 scalar_r = real(scalar_d, kind=rp)
382
383 call signed_distance_field(temp_field, boundary_mesh, scalar_d)
384 call step_function_field(temp_field, scalar_r, 1.0_rp, 0.0_rp)
385
386 case default
387 call neko_error('Unknown distance transform')
388 end select
389
390 ! Update the global indicator field by max operator
391 if (neko_bcknd_device .eq. 1) then
392 call device_pwmax(this%indicator%x_d, temp_field%x_d, &
393 this%indicator%size())
394 else
395 this%indicator%x = max(this%indicator%x, temp_field%x)
396 end if
397
398 end subroutine init_boundary_mesh
399
401 subroutine init_point_zone(this, json)
402 class(brinkman_source_term_t), intent(inout) :: this
403 type(json_file), intent(inout) :: json
404
405 ! Options
406 character(len=:), allocatable :: zone_name
407
408 type(field_t) :: temp_field
409 class(point_zone_t), pointer :: my_point_zone
410 integer :: i
411
412 ! ------------------------------------------------------------------------ !
413 ! Read the options for the point zone
414
415 call json_get(json, 'name', zone_name)
416
417 ! Compute the indicator field
418 call temp_field%init(this%coef%dof)
419
420 my_point_zone => neko_point_zone_registry%get_point_zone(zone_name)
421
422 if (neko_bcknd_device .eq. 1) then
423 call device_cfill_mask(temp_field%x_d, 1.0_rp, temp_field%size(), &
424 my_point_zone%mask_d, my_point_zone%size)
425 else
426 call cfill_mask(temp_field%x, 1.0_rp, temp_field%size(), &
427 my_point_zone%mask, my_point_zone%size)
428 end if
429
430 ! Update the global indicator field by max operator
431 if (neko_bcknd_device .eq. 1) then
432 call device_pwmax(this%indicator%x_d, temp_field%x_d, &
433 this%indicator%size())
434 else
435 this%indicator%x = max(this%indicator%x, temp_field%x)
436 end if
437
438 end subroutine init_point_zone
439
440end module brinkman_source_term
double real
Copy data between host and device (or device and device)
Definition device.F90:65
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_compute(this, t, tstep)
Computes the source term and adds the result to fields.
subroutine brinkman_source_term_free(this)
Destructor.
subroutine brinkman_source_term_init_from_json(this, json, fields, coef)
The common constructor using a JSON object.
Coefficients.
Definition coef.f90:34
subroutine, public device_cfill_mask(a_d, c, size, mask_d, mask_size)
Fill a constant to a masked vector. .
Device abstraction, common interface for various accelerators.
Definition device.F90:34
integer, parameter, public host_to_device
Definition device.F90:46
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
A module containing filter functions and subroutines. These functions are used to modify fields in a ...
Definition filters.f90:36
subroutine, public permeability_field(f, k_0, k_1, q)
Apply a permeability function to a field.
Definition filters.f90:90
subroutine, public smooth_step_field(f, edge0, edge1)
Apply a smooth step function to a field.
Definition filters.f90:70
subroutine, public step_function_field(f, x0, value0, value1)
Apply a step function to a field.
Definition filters.f90:109
Implements fld_file_output_t.
Utilities for retrieving parameters from the case files.
Logging routines.
Definition log.f90:34
type(log_t), public neko_log
Global log stream.
Definition log.f90:65
integer, parameter, public log_size
Definition log.f90:42
Definition math.f90:60
subroutine, public cfill_mask(a, c, size, mask, mask_size)
Fill a constant to a masked vector. .
Definition math.f90:296
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 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:54
Base abstract class for filter.
Definition filter.f90:48
A simple output saving a list of fields to a .fld file.
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.
#define max(a, b)
Definition tensor.cu:40