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