Neko  0.9.0
A portable framework for high-order spectral element flow simulations
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 num_types, only: rp, dp
36  use field, only: field_t
37  use field_list, only: field_list_t
38  use json_module, only: json_file
41  use source_term, only: source_term_t
42  use coefs, only: coef_t
44  use utils, only: neko_error
45  use field_math, only: field_subcol3
46 
47  use math, only: pwmax
48  use device_math, only: device_pwmax
49  implicit none
50  private
51 
55  type, public, extends(source_term_t) :: brinkman_source_term_t
56  private
57 
59  type(field_t) :: indicator
61  type(field_t) :: brinkman
62  contains
64  procedure, public, pass(this) :: init => &
67  procedure, public, pass(this) :: free => brinkman_source_term_free
69  procedure, public, pass(this) :: compute_ => brinkman_source_term_compute
70 
71  ! ----------------------------------------------------------------------- !
72  ! Private methods
73  procedure, pass(this) :: init_boundary_mesh
74  procedure, pass(this) :: init_point_zone
75  end type brinkman_source_term_t
76 
77 contains
78 
79  ! ========================================================================== !
80  ! Public methods
81 
86  subroutine brinkman_source_term_init_from_json(this, json, fields, coef)
87  use file, only: file_t
88  use tri_mesh, only: tri_mesh_t
94  use json_module, only: json_core, json_value
95  implicit none
96 
97  class(brinkman_source_term_t), intent(inout) :: this
98  type(json_file), intent(inout) :: json
99  type(field_list_t), intent(inout), target :: fields
100  type(coef_t), intent(inout), target :: coef
101  real(kind=rp) :: start_time, end_time
102 
103  character(len=:), allocatable :: filter_type
104  real(kind=rp), dimension(:), allocatable :: brinkman_limits
105  real(kind=rp) :: brinkman_penalty
106 
107  type(json_value), pointer :: json_object_list
108  type(json_core) :: core
109 
110  character(len=:), allocatable :: object_type
111  type(json_file) :: object_settings
112  integer :: n_regions
113  integer :: i
114 
115  ! Mandatory fields for the general source term
116  call json_get_or_default(json, "start_time", start_time, 0.0_rp)
117  call json_get_or_default(json, "end_time", end_time, huge(0.0_rp))
118 
119  ! Read the options for the permeability field
120  call json_get(json, 'brinkman.limits', brinkman_limits)
121  call json_get(json, 'brinkman.penalty', brinkman_penalty)
122 
123  if (size(brinkman_limits) .ne. 2) then
124  call neko_error('brinkman_limits must be a 2 element array of reals')
125  end if
126 
127  call this%free()
128  call this%init_base(fields, coef, start_time, end_time)
129 
130  ! ------------------------------------------------------------------------ !
131  ! Allocate the permeability and indicator field
132 
133  if (neko_field_registry%field_exists('brinkman_indicator') &
134  .or. neko_field_registry%field_exists('brinkman')) then
135  call neko_error('Brinkman field already exists.')
136  end if
137 
138  call this%indicator%init(coef%dof)
139  call this%brinkman%init(coef%dof)
140 
141  ! ------------------------------------------------------------------------ !
142  ! Select which constructor should be called
143 
144  call json%get('objects', json_object_list)
145  call json%info('objects', n_children = n_regions)
146  call json%get_core(core)
147 
148  do i = 1, n_regions
149  call json_extract_item(core, json_object_list, i, object_settings)
150  call json_get_or_default(object_settings, 'type', object_type, 'none')
151 
152  select case (object_type)
153  case ('boundary_mesh')
154  call this%init_boundary_mesh(object_settings)
155  case ('point_zone')
156  call this%init_point_zone(object_settings)
157 
158  case ('none')
159  call object_settings%print()
160  call neko_error('Brinkman source term objects require a region type')
161  case default
162  call neko_error('Brinkman source term unknown region type')
163  end select
164 
165  end do
166 
167  ! Run filter on the full indicator field to smooth it out.
168  call json_get_or_default(json, 'filter.type', filter_type, 'none')
169 
170  select case (filter_type)
171  case ('none')
172  ! Do nothing
173  case default
174  call neko_error('Brinkman source term unknown filter type')
175  end select
176 
177  ! ------------------------------------------------------------------------ !
178  ! Compute the permeability field
179 
180  this%brinkman = this%indicator
181  call permeability_field(this%brinkman, &
182  brinkman_limits(1), brinkman_limits(2), brinkman_penalty)
183 
185 
187  subroutine brinkman_source_term_free(this)
188  class(brinkman_source_term_t), intent(inout) :: this
189 
190  call this%indicator%free()
191  call this%brinkman%free()
192  call this%free_base()
193  end subroutine brinkman_source_term_free
194 
198  subroutine brinkman_source_term_compute(this, t, tstep)
199  class(brinkman_source_term_t), intent(inout) :: this
200  real(kind=rp), intent(in) :: t
201  integer, intent(in) :: tstep
202  type(field_t), pointer :: u, v, w, fu, fv, fw
203  integer :: n
204 
205  n = this%fields%item_size(1)
206 
207  u => neko_field_registry%get_field('u')
208  v => neko_field_registry%get_field('v')
209  w => neko_field_registry%get_field('w')
210 
211  fu => this%fields%get(1)
212  fv => this%fields%get(2)
213  fw => this%fields%get(3)
214 
215  call field_subcol3(fu, u, this%brinkman, n)
216  call field_subcol3(fv, v, this%brinkman, n)
217  call field_subcol3(fw, w, this%brinkman, n)
218 
219  end subroutine brinkman_source_term_compute
220 
221  ! ========================================================================== !
222  ! Private methods
223 
225  subroutine init_boundary_mesh(this, json)
226  use file, only: file_t
227  use tri_mesh, only: tri_mesh_t
233  use aabb, only : aabb_t, get_aabb
234  implicit none
235 
236  class(brinkman_source_term_t), intent(inout) :: this
237  type(json_file), intent(inout) :: json
238 
239  ! Options
240  character(len=:), allocatable :: mesh_file_name
241  character(len=:), allocatable :: distance_transform
242  character(len=:), allocatable :: filter_type
243  character(len=:), allocatable :: mesh_transform
244 
245  ! Read the options for the boundary mesh
246  type(file_t) :: mesh_file
247  type(tri_mesh_t) :: boundary_mesh
248  real(kind=rp) :: scalar_r
249  real(kind=dp) :: scalar_d
250 
251  ! Mesh transform options variables
252  real(kind=dp), dimension(:), allocatable :: box_min, box_max
253  logical :: keep_aspect_ratio
254  real(kind=dp), dimension(3) :: scaling
255  real(kind=dp), dimension(3) :: translation
256  type(field_t) :: temp_field
257  type(aabb_t) :: mesh_box, target_box
258  integer :: idx_p
259 
260  ! ------------------------------------------------------------------------ !
261  ! Read the options for the boundary mesh
262 
263  call json_get(json, 'name', mesh_file_name)
264 
265  ! Settings on how to filter the design field
266  call json_get(json, 'distance_transform.type', distance_transform)
267 
268  ! ------------------------------------------------------------------------ !
269  ! Load the immersed boundary mesh
270 
271  mesh_file = file_t(mesh_file_name)
272  call mesh_file%read(boundary_mesh)
273 
274  if (boundary_mesh%nelv .eq. 0) then
275  call neko_error('No elements in the boundary mesh')
276  end if
277 
278  ! ------------------------------------------------------------------------ !
279  ! Transform the mesh if specified.
280 
281  call json_get_or_default(json, 'mesh_transform.type', &
282  mesh_transform, 'none')
283 
284  select case (mesh_transform)
285  case ('none')
286  ! Do nothing
287  case ('bounding_box')
288  call json_get(json, 'mesh_transform.box_min', box_min)
289  call json_get(json, 'mesh_transform.box_max', box_max)
290  call json_get_or_default(json, 'mesh_transform.keep_aspect_ratio', &
291  keep_aspect_ratio, .true.)
292 
293  if (size(box_min) .ne. 3 .or. size(box_max) .ne. 3) then
294  call neko_error('Case file: mesh_transform. &
295  &box_min and box_max must be 3 element arrays of reals')
296  end if
297 
298  call target_box%init(box_min, box_max)
299 
300  mesh_box = get_aabb(boundary_mesh)
301 
302  scaling = target_box%get_diagonal() / mesh_box%get_diagonal()
303  if (keep_aspect_ratio) then
304  scaling = minval(scaling)
305  end if
306 
307  translation = - scaling * mesh_box%get_min() + target_box%get_min()
308 
309  do idx_p = 1, boundary_mesh%mpts
310  boundary_mesh%points(idx_p)%x = &
311  scaling * boundary_mesh%points(idx_p)%x + translation
312  end do
313 
314  case default
315  call neko_error('Unknown mesh transform')
316  end select
317 
318  ! ------------------------------------------------------------------------ !
319  ! Compute the permeability field
320 
321  ! Assign the signed distance field to all GLL points in the permeability
322  ! field. Initally we just run a brute force loop over all GLL points and
323  ! compute the signed distance function. This should be replaced with a
324  ! more efficient method, such as a tree search.
325 
326  call temp_field%init(this%coef%dof)
327 
328  ! Select how to transform the distance field to a design field
329  select case (distance_transform)
330  case ('smooth_step')
331  call json_get(json, 'distance_transform.value', scalar_d)
332  scalar_r = real(scalar_d, kind=rp)
333 
334  call signed_distance_field(temp_field, boundary_mesh, scalar_d)
335  call smooth_step_field(temp_field, scalar_r, 0.0_rp)
336 
337  case ('step')
338 
339  call json_get(json, 'distance_transform.value', scalar_d)
340  scalar_r = real(scalar_d, kind=rp)
341 
342  call signed_distance_field(temp_field, boundary_mesh, scalar_d)
343  call step_function_field(temp_field, scalar_r, 1.0_rp, 0.0_rp)
344 
345  case default
346  call neko_error('Unknown distance transform')
347  end select
348 
349  ! ------------------------------------------------------------------------ !
350  ! Run filter on the temporary indicator field to smooth it out.
351  call json_get_or_default(json, 'filter.type', filter_type, 'none')
352 
353  select case (filter_type)
354  case ('none')
355  ! Do nothing
356  case default
357  call neko_error('Unknown filter type')
358  end select
359 
360  ! Update the global indicator field by max operator
361  if (neko_bcknd_device .eq. 1) then
362  call device_pwmax(this%indicator%x_d, temp_field%x_d, &
363  this%indicator%size())
364  else
365  this%indicator%x = max(this%indicator%x, temp_field%x)
366  end if
367 
368  end subroutine init_boundary_mesh
369 
371  subroutine init_point_zone(this, json)
376  use point_zone, only: point_zone_t
379  implicit none
380 
381  class(brinkman_source_term_t), intent(inout) :: this
382  type(json_file), intent(inout) :: json
383 
384  ! Options
385  character(len=:), allocatable :: zone_name
386  character(len=:), allocatable :: filter_type
387 
388  type(field_t) :: temp_field
389  class(point_zone_t), pointer :: my_point_zone
390  integer :: i
391 
392  ! ------------------------------------------------------------------------ !
393  ! Read the options for the point zone
394 
395  call json_get(json, 'name', zone_name)
396  call json_get_or_default(json, 'filter.type', filter_type, 'none')
397 
398  ! Compute the indicator field
399  call temp_field%init(this%coef%dof)
400 
401  my_point_zone => neko_point_zone_registry%get_point_zone(zone_name)
402 
403  do i = 1, my_point_zone%size
404  temp_field%x(my_point_zone%mask(i), 1, 1, 1) = 1.0_rp
405  end do
406 
407  ! Run filter on the temporary indicator field to smooth it out.
408 
409  select case (filter_type)
410  case ('none')
411  ! Do nothing
412  case default
413  call neko_error('Unknown filter type')
414  end select
415 
416  ! Update the global indicator field by max operator
417  this%indicator%x = max(this%indicator%x, temp_field%x)
418 
419  end subroutine init_point_zone
420 
421 end module brinkman_source_term
double real
Definition: device_config.h:12
Copy data between host and device (or device and device)
Definition: device.F90:51
Retrieves a parameter by name or assigns a provided default value. In the latter case also adds the m...
Definition: json_utils.f90:54
Retrieves a parameter by name or throws an error.
Definition: json_utils.f90:45
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
Device abstraction, common interface for various accelerators.
Definition: device.F90:34
integer, parameter, public host_to_device
Definition: device.F90:47
subroutine, public field_subcol3(a, b, c, n)
Returns .
Definition: field_math.f90:517
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
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
Utilities for retrieving parameters from the case files.
Definition: json_utils.f90:34
Definition: math.f90:60
Build configurations.
Definition: neko_config.f90:34
integer, parameter neko_bcknd_device
Definition: neko_config.f90:44
integer, parameter, public dp
Definition: num_types.f90:9
integer, parameter, public rp
Global precision used in computations.
Definition: num_types.f90:12
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.
Definition: source_term.f90:34
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
Definition: field_list.f90:13
A wrapper around a polymorphic generic_file_t that handles its init. This is essentially a factory fo...
Definition: file.f90:54
Base abstract type for point zones.
Definition: point_zone.f90:47
Base abstract type for source terms.
Definition: source_term.f90:43
#define max(a, b)
Definition: tensor.cu:40