Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
design_brinkman.f90
1! Copyright (c) 2024, 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
33! Implements the `brinkman_design_t` type.
34module brinkman_design
35 use num_types, only: rp, sp
36 use field, only: field_t
37 use json_module, only: json_file
39 use coefs, only: coef_t
41 use scratch_registry, only: neko_scratch_registry
42 use fld_file_output, only: fld_file_output_t
43 use point_zone_registry, only: neko_point_zone_registry
44 use point_zone, only: point_zone_t
46 use neko_config, only: neko_bcknd_device
47 use device, only: device_memcpy, host_to_device
48 use design, only: design_t
49 use math, only: rzero
50 use simulation_m, only: simulation_t
51 use json_module, only: json_file
53 use vector, only: vector_t
54 use math, only: copy
55 use device_math, only: device_copy
56 use field_registry, only: neko_field_registry
59 use field_math, only: field_rzero
60 use json_utils, only: json_get, json_get_or_default, json_extract_object
61 use utils, only: neko_error
62 implicit none
63 private
64
66 type, extends(design_t), public :: brinkman_design_t
67 private
68
69 ! TODO
70 ! in the future make this a derived type of a `design_variable`
71 ! type, public, extends(design_variable_t) :: brinkman_design_t
73 type(field_t), pointer :: design_indicator
74
76 ! TODO
77 ! NOTE: Tim, right now we only have the brinkman term
78 ! in the future we may map to other coeeficients for other equations...
79 ! in which case this should be a field list
80 !
81 ! or as I describe below, we also have multiple constraints,
82 ! so a list-of-lists may be the correct way forward
83 type(field_t), pointer :: brinkman_amplitude
84
85 ! NOTE:
86 ! again, we have to be so clear with nomenclature.
87 ! If we have an objective function F.
88 ! I also like to use \rho to describe the design_indicator
89 !
90 ! and let's say, for completness, we're doing conjugate heat transfer,
91 ! so we have to map
92 ! to 3 coefficients, \chi, C and \kappa.
93 !
94 ! then we will perform 3 mapping,
95 ! \rho -> \chi
96 ! \rho -> C
97 ! \rho -> \kappa
98 !
99 ! Then during the forward/adjoint looping there will be an additional
100 ! object, the "objective_function" object that will be responsible for
101 ! computing the sensitivity of the objective function with respect to the
102 ! coefficients
103 ! ie,
104 ! dF/d\chi, dF/dC and dF/d\kappa
105 !
106 ! What I'm calling "sensitivity" here, is the sensitivity with respect to
107 ! the design indicator
108 ! so dF/d\rho
109 !
110 ! so the proceedure "map_backwards" will take in the field list
111 ! dF/d\chi, dF/dC and dF/d\kappa
112 !and chain rule it's way back to
113 ! dF/d\rho
114 ! and store it here v
115 type(field_t), pointer :: sensitivity
116 ! have a field list here
117 ! type(filed_list_t), public :: constraint_sensitivity
118 ! HOWEVER !
119 ! What is a bit confusing to me... is how we'll deal with constraints.
120 !
121 ! Because in principle we could have constraints C1, C2, C3 etc
122 ! Implying we also want dC1/d\rho, dC2/d\rho etc
123 ! So this "sensitivity" should also be a list...
124 !
125 ! So I bet you'll have a nice abstract way of constructing this in the
126 ! future but I think a sort of list-of-lists would be nice.
127 !
128 ! For F:
129 ! \rho -> \tild{\rho} -> \chi
130 ! \rho -> \tild{\rho} -> C
131 ! \rho -> \tild{\rho} -> \kappa
132 !
133 ! For C1:
134 ! \rho -> \tild{\rho} (eg for a volume constraint or so)
135 !
136 ! For C2:
137 ! \rho -> \tild{\rho} (for something else)
138 !
139 ! etc..
140 !
141 ! So now we have multiple instances of the "objective" type,
142 ! each for F, C1, C2 etc
143 ! each of them can pass their dF\d_coefficents to the design
144 ! and we continue from there.
145 !
146 ! perhaps even the "objective" type is defined as a list of objectives.
147
148 ! Let's say way have a chain of two mappings
149
150 ! This is still up for discussion where the mapping will eventually live.
151 ! Maybe it makes more sense to keep the design as clean as possible,
152 ! and then the resposibility of converting the design into coefficients
153 ! used in the PDE being solved would be the responsibility of the module
154 ! used for solving the PDE, ie, the simulation.
155 ! Then the simulation would responsiblity for asking "do we have a scalar?"
156 ! ok, let me figure out how to map this design into a coeffient in the
157 ! scalar equation.
158 !
159 ! I guess then the mapping backwards gets confusing to me.
160 ! anyway, for now, let's keep things how they are and I suppose the
161 ! intention is to have different types of designs for different types of
162 ! problems.
163
169
171 class(point_zone_t), pointer :: optimization_domain
173 logical :: has_mask
174
175 ! TODO
176 ! you also had logicals for convergence etc,
177 ! I feel they should live in "problem"
178
179 ! afield writer would be nice too
180 type(fld_file_output_t), private :: output
181
182 contains
183
184 ! ----------------------------------------------------------------------- !
185 ! Initializations
186
188 generic, public :: init => init_from_json_sim, init_from_components
190 procedure, pass(this), public :: init_from_json_sim => &
191 brinkman_design_init_from_json_sim
193 procedure, pass(this), public :: init_from_components => &
194 brinkman_design_init_from_components
195
197 procedure, pass(this) :: get_values => brinkman_design_get_design
198
200 procedure, pass(this) :: design_get_x => brinkman_design_get_x
202 procedure, pass(this) :: design_get_x_i => brinkman_design_get_x_i
204 procedure, pass(this) :: design_get_y => brinkman_design_get_y
206 procedure, pass(this) :: design_get_y_i => brinkman_design_get_y_i
208 procedure, pass(this) :: design_get_z => brinkman_design_get_z
210 procedure, pass(this) :: design_get_z_i => brinkman_design_get_z_i
211
213 procedure, pass(this) :: update_design => brinkman_design_update_design
214
216 ! design_indicator -> filtering -> chi
217 ! and ultimately handle mapping different coeficients!
218 procedure, pass(this) :: map_forward => brinkman_design_map_forward
220 ! d_design_indicator <- d_filtering <- d_chi
221 ! and ultimately handle mapping different coeficients!
222 procedure, pass(this) :: map_backward => brinkman_design_map_backward
223 ! TODO
224 ! maybe it would have been smarter to have a "coeficient" type,
225 ! which is just a scalar field and set of mappings going from
226 ! design_indicator -> coeficient and their corresponding chain rules
227 ! maybe also some information about what equation they live in...
228
229 ! a writer being called from outside would be nice
230 procedure, pass(this) :: write => brinkman_design_write
231
233 procedure, pass(this) :: free => brinkman_design_free
234 ! TODO
235 ! I'm not sure who owns the optimizer...
236 ! but it would make sense to have it in here so you provide it
237 ! with dF/d_design_indicator and it updates itself.
238 ! procedure, pass(this) :: update => brinkman_design_update_design
239 end type brinkman_design_t
240
241contains
242
244 subroutine brinkman_design_init_from_json_sim(this, parameters, simulation)
245 class(brinkman_design_t), intent(inout) :: this
246 type(json_file), intent(inout) :: parameters
247 type(simulation_t), intent(inout) :: simulation
248 type(json_file) :: json_subdict
249 character(len=:), allocatable :: domain_name, domain_type, name
250
251 call json_get_or_default(parameters, 'name', name, 'Brinkman Design')
252 call json_get_or_default(parameters, 'domain.type', domain_type, 'full')
253
254 select case (trim(domain_type))
255 case ('full')
256 this%has_mask = .false.
257 case ('point_zone')
258 this%has_mask = .true.
259 call json_get(parameters, 'domain.zone_name', domain_name)
260 this%optimization_domain => &
261 neko_point_zone_registry%get_point_zone(domain_name)
262
263 case default
264 call neko_error('brinkman design only supports point_zones for ' // &
265 'optimization domain types')
266
267 end select
268
269 ! Initialize and inject into the simulation
270 call this%init_from_components(name, simulation)
271
272 ! Initialize the mapper
273 associate(coef => simulation%neko_case%fluid%c_Xh, &
274 gs => simulation%neko_case%fluid%gs_Xh)
275
276 if ('mapping' .in. parameters) then
277 call json_extract_object(parameters, 'mapping', json_subdict)
278 call this%mapping%init_base(coef)
279 call this%mapping%add(parameters, 'mapping')
280 end if
281
282 if ('initial_distribution' .in. parameters) then
283 call json_extract_object(parameters, 'initial_distribution', json_subdict)
284 call set_optimization_ic(this%design_indicator, coef, gs, &
285 json_subdict)
286 else
287 call field_rzero(this%design_indicator)
288 end if
289 end associate
290
291 ! Map to the Brinkman amplitude
292 call this%map_forward()
293
294 end subroutine brinkman_design_init_from_json_sim
295
297 subroutine brinkman_design_free(this)
298 class(brinkman_design_t), intent(inout) :: this
299
300 call this%free_base()
301 call this%brinkman_amplitude%free()
302 call this%design_indicator%free()
303 call this%sensitivity%free()
304
305 end subroutine brinkman_design_free
306
307 subroutine brinkman_design_init_from_components(this, name, simulation)
308 class(brinkman_design_t), intent(inout) :: this
309 character(len=*), intent(in) :: name
310 type(simulation_t), intent(inout) :: simulation
311 integer :: n, i
312 type(simple_brinkman_source_term_t) :: forward_brinkman, adjoint_brinkman
313
314 associate(dof => simulation%neko_case%fluid%dm_Xh)
315
316 call neko_field_registry%add_field(dof, "design_indicator", .true.)
317 call neko_field_registry%add_field(dof, "brinkman_amplitude", .true.)
318 call neko_field_registry%add_field(dof, "sensitivity", .true.)
319
320 end associate
321
322 this%design_indicator => &
323 neko_field_registry%get_field("design_indicator")
324 this%brinkman_amplitude => &
325 neko_field_registry%get_field("brinkman_amplitude")
326 this%sensitivity => &
327 neko_field_registry%get_field("sensitivity")
328
329 ! TODO
330 ! this is where we steal basically everything in
331 ! brinkman_source_term regarding loading initial fields
332 ! for now, make it a cylinder by hand
333 this%design_indicator = 0.0_rp
334 this%brinkman_amplitude = 0.0_rp
335 this%design_indicator%x = 0.0_rp
336
337 n = this%design_indicator%dof%size()
338 ! This is probably getting fixed in tim's PR anyway, otherwise I'll fix it.
339 do i = 1, n
340 this%design_indicator%x(i,1,1,1) = 0.0_rp
341 end do
342
343 ! again this will be handled better in the future...
344 if (neko_bcknd_device .eq. 1) then
345 call device_memcpy(this%design_indicator%x, &
346 this%design_indicator%x_d, n, &
347 host_to_device, sync = .false.)
348 end if
349
350 ! TODO
351 ! Regarding masks and filters,
352 ! I suppose there are two ways of thinking about it:
353 ! 1) Mask first, then filter
354 ! 2) filter first, then mask
355 !
356 ! Each one can be used to achieve different results, and when we do complex
357 ! non-linear filter cascades the choice here has implications for exactly
358 ! how we define "minimum size control".
359 !
360 ! On top of this, I've only ever used spatial convolution filters before,
361 ! not PDE based filters.
362 ! With these filters you need to think of how your domain is "padded" when
363 ! you move the kernel over the boundary. Again, you can achieve different
364 ! effects depending on the choice of padding.
365 ! I don't really know if there's a similar implication for PDE based
366 ! based filters, I bet there is. (We should ask Niels and Casper) I bet it
367 ! means we need to consider the boundary of the mask as the boundary we
368 ! enforce to be Nuemann, not the boundary of the computational domain.
369 ! That sounds REALLY hard to implement...I hope it doesn't come down to
370 ! that.
371 !
372 ! We can always change this decision later, but I'm going with (1)
373 ! mask first, then filter.
374 ! The reason being, one of the purposes of the filtering is to avoid sharp
375 ! interfaces, if we filter first then mask there's a chance we have a sharp
376 ! interface on the boundary of the optimization domain.
377 ! if we mask first then filter, at least all the boundaries will be smooth.
378 if (this%has_mask) then
379 call mask_exterior_const(this%design_indicator, &
380 this%optimization_domain, 0.0_rp)
381 end if
382
383 ! a field writer would be nice to output
384 ! - design indicator (\rho)
385 ! - mapped design (\chi)
386 ! - sensitivity (dF/d\chi)
387 ! TODO
388 ! do this properly with JSON
389 ! TODO
390 ! obviously when we do the mappings properly, to many coeficients,
391 ! we'll also have to modify this
392 call this%output%init(sp, 'design', 3)
393 call this%output%fields%assign_to_field(1, this%design_indicator)
394 call this%output%fields%assign_to_field(2, this%brinkman_amplitude)
395 call this%output%fields%assign_to_field(3, this%sensitivity)
396
397 call this%init_base(name, n)
398
399 ! init the simple brinkman term for the forward problem
400 call forward_brinkman%init_from_components( &
401 simulation%fluid%f_x, &
402 simulation%fluid%f_y, &
403 simulation%fluid%f_z, &
404 this%brinkman_amplitude, &
405 simulation%fluid%u, &
406 simulation%fluid%v, &
407 simulation%fluid%w, &
408 simulation%fluid%c_Xh)
409 ! append brinkman source term to the forward problem
410 call simulation%fluid%source_term%add(forward_brinkman)
411
412 ! init the simple brinkman term for the adjoint
413 call adjoint_brinkman%init_from_components( &
414 simulation%adjoint_fluid%f_adj_x, &
415 simulation%adjoint_fluid%f_adj_y, &
416 simulation%adjoint_fluid%f_adj_z, &
417 this%brinkman_amplitude, &
418 simulation%adjoint_fluid%u_adj, &
419 simulation%adjoint_fluid%v_adj, &
420 simulation%adjoint_fluid%w_adj, &
421 simulation%adjoint_fluid%c_Xh)
422 ! append brinkman source term based on design
423
424 select type (f => simulation%adjoint_fluid)
425 type is (adjoint_fluid_pnpn_t)
426 call f%source_term%add(adjoint_brinkman)
427 class default
428 end select
429
430 end subroutine brinkman_design_init_from_components
431
432
433 subroutine brinkman_design_map_forward(this)
434 class(brinkman_design_t), intent(inout) :: this
435
436 ! TODO, see previous todo about mask first, then mapping
437 if (this%has_mask) then
438 call mask_exterior_const(this%design_indicator, &
439 this%optimization_domain, 0.0_rp)
440 end if
441
442 call this%mapping%apply_forward(this%brinkman_amplitude, &
443 this%design_indicator)
444
445 end subroutine brinkman_design_map_forward
446
447 subroutine brinkman_design_get_design(this, values)
448 class(brinkman_design_t), intent(in) :: this
449 type(vector_t), intent(inout) :: values
450 integer :: n
451
452 n = this%size()
453 call values%init(n)
454 call copy(values%x, this%design_indicator%x, n)
455 if (neko_bcknd_device .eq. 1) then
456 call device_copy(values%x_d, this%design_indicator%x_d, n)
457 end if
458
459 end subroutine brinkman_design_get_design
460
461 subroutine brinkman_design_get_x(this, x)
462 class(brinkman_design_t), intent(in) :: this
463 type(vector_t), intent(inout) :: x
464 integer :: n
465
466 n = this%size()
467 call x%init(n)
468 call copy(x%x, this%design_indicator%dof%x, n)
469 if (neko_bcknd_device .eq. 1) then
470 call device_copy(x%x_d, this%design_indicator%dof%x_d, n)
471 end if
472
473 end subroutine brinkman_design_get_x
474
475 function brinkman_design_get_x_i(this, i) result(x_i)
476 class(brinkman_design_t), intent(in) :: this
477 integer, intent(in) :: i
478 real(kind=rp) :: x_i
479 integer :: n
480
481 n = this%size()
482 if (i .lt. 1 .or. i .gt. n) then
483 call neko_error('brinkman_design_get_x_i: index out of bounds')
484 end if
485
486 x_i = this%design_indicator%dof%x(i,1,1,1)
487
488 end function brinkman_design_get_x_i
489
490 subroutine brinkman_design_get_y(this, y)
491 class(brinkman_design_t), intent(in) :: this
492 type(vector_t), intent(inout) :: y
493 integer :: n
494
495 n = this%size()
496 call y%init(n)
497 call copy(y%x, this%design_indicator%dof%y, n)
498 if (neko_bcknd_device .eq. 1) then
499 call device_copy(y%x_d, this%design_indicator%dof%y_d, n)
500 end if
501
502 end subroutine brinkman_design_get_y
503
504 function brinkman_design_get_y_i(this, i) result(y_i)
505 class(brinkman_design_t), intent(in) :: this
506 integer, intent(in) :: i
507 real(kind=rp) :: y_i
508 integer :: n
509
510 n = this%size()
511 if (i .lt. 1 .or. i .gt. n) then
512 call neko_error('brinkman_design_get_y_i: index out of bounds')
513 end if
514
515 y_i = this%design_indicator%dof%y(i,1,1,1)
516
517 end function brinkman_design_get_y_i
518
519 subroutine brinkman_design_get_z(this, z)
520 class(brinkman_design_t), intent(in) :: this
521 type(vector_t), intent(inout) :: z
522 integer :: n
523
524 n = this%size()
525 call z%init(n)
526 call copy(z%x, this%design_indicator%dof%z, n)
527 if (neko_bcknd_device .eq. 1) then
528 call device_copy(z%x_d, this%design_indicator%dof%z_d, n)
529 end if
530
531 end subroutine brinkman_design_get_z
532
533 function brinkman_design_get_z_i(this, i) result(z_i)
534 class(brinkman_design_t), intent(in) :: this
535 integer, intent(in) :: i
536 real(kind=rp) :: z_i
537 integer :: n
538
539 n = this%size()
540 if (i .lt. 1 .or. i .gt. n) then
541 call neko_error('brinkman_design_get_z_i: index out of bounds')
542 end if
543
544 z_i = this%design_indicator%dof%z(i,1,1,1)
545
546 end function brinkman_design_get_z_i
547
548 subroutine brinkman_design_update_design(this, values)
549 class(brinkman_design_t), intent(inout) :: this
550 type(vector_t), intent(inout) :: values
551 integer :: n
552
553 n = this%size()
554 call copy(this%design_indicator%x, values%x, n)
555 if (neko_bcknd_device .eq. 1) then
556 call device_copy(this%design_indicator%x_d, values%x_d, n)
557 end if
558
559 call this%map_forward()
560
561 call copy(values%x, this%design_indicator%x, n)
562 if (neko_bcknd_device .eq. 1) then
563 call device_copy(values%x_d, this%design_indicator%x_d, n)
564 end if
565
566 end subroutine brinkman_design_update_design
567
568 subroutine brinkman_design_map_backward(this, sensitivity)
569 class(brinkman_design_t), intent(inout) :: this
570 type(vector_t), intent(in) :: sensitivity
571 type(field_t), pointer :: tmp_fld
572 integer :: temp_indices(1)
573
574 call neko_scratch_registry%request_field(tmp_fld, temp_indices(1))
575
576 call vector_to_field(tmp_fld, sensitivity)
577
578 call this%mapping%apply_backward(this%sensitivity, tmp_fld)
579
580 ! TODO
581 ! DELETE THIS LATER
582 !
583 ! When Abbas writes the interface for the optimization
584 ! module this may be a moot point, because we would only really collect
585 ! the sensitivity of the design variables inside the mask.
586 !
587 ! Note for Abbas,
588 ! I'm NOT doing this because I'm too lazy and I just need masks so I can
589 ! test something in the passive scalar.
590 if (this%has_mask) then
591 call mask_exterior_const(this%sensitivity, this%optimization_domain, &
592 0.0_rp)
593 end if
594
595 call neko_scratch_registry%relinquish_field(temp_indices)
596
597 end subroutine brinkman_design_map_backward
598
599 subroutine brinkman_design_write(this, idx)
600 class(brinkman_design_t), intent(inout) :: this
601 integer, intent(in) :: idx
602
603 call this%output%sample(real(idx, kind=rp))
604
605 end subroutine brinkman_design_write
606
607end module brinkman_design
Adjoint Pn/Pn formulation.
Implements the design_t.
Definition design.f90:34
Implements the mapping_handler_t type.
Mappings to be applied to a scalar field.
Definition mapping.f90:35
Some common Masking operations we may need.
Definition mask_ops.f90:34
Contains extensions to the neko library required to run the topology optimization code.
Definition neko_ext.f90:9
subroutine, public field_to_vector(vector, field)
Field to vector.
Definition neko_ext.f90:256
subroutine, public vector_to_field(field, vector)
Vector to field.
Definition neko_ext.f90:232
Optimization initial condition.
Implements the simple_brinkman_source_term_t type.
Implements the steady_problem_t type.
A topology optimization design variable.
An abstract design type.
Definition design.f90:52
Abstract class for handling mapping_cascade.