Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
set_optimization_ic.f90
1! Copyright (c) 2021, 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 gather_scatter, only : gs_t, gs_op_add
36 use neko_config, only : neko_bcknd_device
37 use num_types, only : rp
38 use device_math, only : device_col2
39 use device, only : device_memcpy, host_to_device
40 use field, only : field_t
41 use utils, only : neko_error, filename_chsuffix, filename_suffix, &
42 neko_warning, neko_fname_len, extract_fld_file_index
43 use coefs, only : coef_t
44 use math, only : col2, cfill, cfill_mask
45 use json_module, only : json_file
46 use json_utils, only: json_get, json_get_or_default
47 use point_zone, only: point_zone_t
48 use point_zone_registry, only: neko_point_zone_registry
49 use field_registry, only: neko_field_registry
50 use logger, only: neko_log, log_size
51 use fld_file_data, only: fld_file_data_t
52 use fld_file, only: fld_file_t
53 use checkpoint, only: chkp_t
54 use file, only: file_t
55 use global_interpolation, only: global_interpolation_t
56 use interpolation, only: interpolator_t
57 use space, only: space_t, gll
58 implicit none
59 private
60
62 module procedure set_optimization_ic_int
63 end interface set_optimization_ic
64
65 public :: set_optimization_ic
66
67contains
68
83 subroutine set_optimization_ic_int(fld, coef, gs, params)
84 type(field_t), intent(inout) :: fld
85 type(coef_t), intent(in) :: coef
86 type(gs_t), intent(inout) :: gs
87 type(json_file), intent(inout) :: params
88
89 ! Variables for retrieving JSON parameters
90 character(len=:), allocatable :: type
91 real(kind=rp) :: ic_value
92 character(len=:), allocatable :: read_str
93 character(len=:), allocatable :: fname, mesh_fname
94 real(kind=rp) :: zone_value, tol
95 logical :: interpolate
96
97 call json_get(params, "type", type)
98
99 select case (trim(type))
100 case ('uniform')
101 call json_get(params, 'value', ic_value)
102 call set_optimization_ic_uniform(fld, ic_value)
103
104 case ('point_zone')
105 call json_get(params, 'base_value', ic_value)
106 call json_get(params, 'zone_name', read_str)
107 call json_get(params, 'zone_value', zone_value)
108
109 call set_optimization_ic_point_zone(fld, ic_value, read_str, zone_value)
110
111 case ('field')
112 call json_get(params, 'file_name', fname)
113 call json_get_or_default(params, 'interpolate', interpolate, .false.)
114 call json_get_or_default(params, 'tolerance', tol, 0.000001_rp)
115 call json_get_or_default(params, 'mesh_file_name', mesh_fname, "none")
116
117 call set_optimization_ic_fld(fld, fname, interpolate, tol, mesh_fname)
118
119 case default
120 call neko_error('Invalid initial condition')
121
122 end select
123
124 call set_optimization_ic_common(fld, coef, gs)
125
126 end subroutine set_optimization_ic_int
127
134 subroutine set_optimization_ic_common(fld, coef, gs)
135 type(field_t), intent(inout) :: fld
136 type(coef_t), intent(in) :: coef
137 type(gs_t), intent(inout) :: gs
138 integer :: n
139
140 n = fld%dof%size()
141 if (neko_bcknd_device .eq. 1) then
142 call device_memcpy(fld%x, fld%x_d, n, &
143 host_to_device, sync = .false.)
144 end if
145
146 ! Ensure continuity across elements for initial conditions
147 call gs%op(fld%x, n, gs_op_add)
148
149 if (neko_bcknd_device .eq. 1) then
150 call device_col2(fld%x_d, coef%mult_d, n)
151 else
152 call col2(fld%x, coef%mult, n)
153 end if
154
155 end subroutine set_optimization_ic_common
156
161 subroutine set_optimization_ic_uniform(fld, ic_value)
162 type(field_t), intent(inout) :: fld
163 real(kind=rp), intent(in) :: ic_value
164 integer :: n
165 character(len=LOG_SIZE) :: log_buf
166
167 call neko_log%message("Type : uniform")
168 write (log_buf, '(A,ES12.6)') "Value: ", ic_value
169 call neko_log%message(log_buf)
170
171 fld = ic_value
172 n = fld%dof%size()
173 if (neko_bcknd_device .eq. 1) then
174 call cfill(fld%x, ic_value, n)
175 end if
176
177 end subroutine set_optimization_ic_uniform
178
186 subroutine set_optimization_ic_point_zone(fld, base_value, zone_name, &
187 zone_value)
188 type(field_t), intent(inout) :: fld
189 real(kind=rp), intent(in) :: base_value
190 character(len=*), intent(in) :: zone_name
191 real(kind=rp), intent(in) :: zone_value
192
193 ! Internal variables
194 character(len=LOG_SIZE) :: log_buf
195 class(point_zone_t), pointer :: zone
196 integer :: size
197
198 call neko_log%message("Type : point_zone")
199 write (log_buf, '(A,ES12.6)') "Base value: ", base_value
200 call neko_log%message(log_buf)
201 call neko_log%message("Zone name : " // trim(zone_name))
202 write (log_buf, '(A,ES12.6)') "Zone value: ", zone_value
203 call neko_log%message(log_buf)
204
205 size = fld%dof%size()
206 zone => neko_point_zone_registry%get_point_zone(trim(zone_name))
207
208 call set_optimization_ic_uniform(fld, base_value)
209 call cfill_mask(fld%x, zone_value, size, zone%mask%get(), zone%size)
210
211 end subroutine set_optimization_ic_point_zone
212
227 subroutine set_optimization_ic_fld(fld, file_name, &
228 interpolate, tolerance, mesh_file_name)
229 type(field_t), intent(inout) :: fld
230 character(len=*), intent(in) :: file_name
231 logical, intent(in) :: interpolate
232 real(kind=rp), intent(in) :: tolerance
233 character(len=*), intent(inout) :: mesh_file_name
234
235 character(len=LOG_SIZE) :: log_buf
236 integer :: sample_idx, sample_mesh_idx
237 type(fld_file_data_t) :: fld_data
238 type(file_t) :: f
239 logical :: mesh_mismatch
240
241 ! ---- For the mesh to mesh interpolation
242 type(global_interpolation_t) :: global_interp
243 ! -----
244
245 ! ---- For space to space interpolation
246 type(space_t) :: prev_Xh
247 type(interpolator_t) :: space_interp
248 ! ----
249
250 call neko_log%message("Type : field")
251 call neko_log%message("File name : " // trim(file_name))
252 write (log_buf, '(A,L1)') "Interpolation: ", interpolate
253 call neko_log%message(log_buf)
254
255 ! Extract sample index from the file name
256 sample_idx = extract_fld_file_index(file_name, -1)
257
258 if (sample_idx .eq. -1) then
259 call neko_error("Invalid file name for the initial condition. The " // &
260 "file format must be e.g. 'mean0.f00001'")
261 end if
262
263 ! Change from "field0.f000*" to "field0.fld" for the fld reader
264 call filename_chsuffix(file_name, trim(file_name), 'fld')
265
266 call fld_data%init()
267 call f%init(trim(file_name))
268
269 if (interpolate) then
270
271 ! If no mesh file is specified, use the default file name
272 if (mesh_file_name .eq. "none") then
273 mesh_file_name = trim(file_name)
274 sample_mesh_idx = sample_idx
275 else
276
277 ! Extract sample index from the mesh file name
278 sample_mesh_idx = extract_fld_file_index(mesh_file_name, -1)
279
280 if (sample_mesh_idx .eq. -1) then
281 call neko_error("Invalid file name for the initial condition." // &
282 " The file format must be e.g. 'mean0.f00001'")
283 end if
284
285 write (log_buf, '(A,ES12.6)') "Tolerance : ", tolerance
286 call neko_log%message(log_buf)
287 write (log_buf, '(A,A)') "Mesh file : ", &
288 trim(mesh_file_name)
289 call neko_log%message(log_buf)
290
291 end if ! if mesh_file_name .eq. none
292
293 ! Read the mesh coordinates if they are not in our fld file
294 if (sample_mesh_idx .ne. sample_idx) then
295 call f%set_counter(sample_mesh_idx)
296 call f%read(fld_data)
297 end if
298
299 end if
300
301 ! Read the field file containing (u,v,w,p)
302 call f%set_counter(sample_idx)
303 call f%read(fld_data)
304
305 !
306 ! Check that the data in the fld file matches the current case.
307 ! Note that this is a safeguard and there are corner cases where
308 ! two different meshes have the same dimension and same # of elements
309 ! but this should be enough to cover obvious cases.
310 !
311 mesh_mismatch = (fld_data%glb_nelv .ne. fld%msh%glb_nelv .or. &
312 fld_data%gdim .ne. fld%msh%gdim)
313
314 if (mesh_mismatch .and. .not. interpolate) then
315 call neko_error("The fld file must match the current mesh! " // &
316 "Use 'interpolate': 'true' to enable interpolation.")
317 else if (.not. mesh_mismatch .and. interpolate) then
318 call neko_log%warning("You have activated interpolation but you " // &
319 "might still be using the same mesh.")
320 end if
321
322
323 ! Mesh interpolation if specified
324 if (interpolate) then
325 ! Issue a warning if the mesh is in single precision
326 select type (ft => f%file_type)
327 type is (fld_file_t)
328 if (.not. ft%dp_precision) then
329 call neko_warning("The coordinates read from the field file " // &
330 "are in single precision.")
331 call neko_log%message("It is recommended to use a mesh in " // &
332 "double precision for better interpolation results.")
333 call neko_log%message("If the interpolation does not work, " // &
334 "you can try to increase the tolerance.")
335 end if
336 class default
337 end select
338
339 ! Generates an interpolator object and performs the point search
340 call fld_data%generate_interpolator(global_interp, fld%dof, fld%msh, &
341 tolerance)
342
343 ! Evaluate design
344 ! NOTE. Is is currently the convention in the design0.f* files that
345 ! the unfiltered design should be stored in the u component.
346 call global_interp%evaluate(fld%x, fld_data%u%x)
347 call global_interp%free()
348
349 else ! No interpolation, just potentially from different spaces
350
351 ! Build a space_t object from the data in the fld file
352 call prev_xh%init(gll, fld_data%lx, fld_data%ly, fld_data%lz)
353 call space_interp%init(fld%Xh, prev_xh)
354
355 ! Do the space-to-space interpolation
356 call space_interp%map_host(fld%x, fld_data%t%x, fld_data%nelv, fld%Xh)
357
358 call space_interp%free()
359
360 end if
361
362 call fld_data%free()
363
364 end subroutine set_optimization_ic_fld
365
366end module optimization_ic
Optimization initial condition.