Neko-TOP
A portable framework for high-order spectral element flow toplogy optimization.
Loading...
Searching...
No Matches
mma_optimizer.f90
1module mma_optimizer
2 use optimizer, only: optimizer_t
3 use problem, only : problem_t
4 use mma, only: mma_t
5 use problem, only: problem_t
6 use num_types, only: rp
7 use utils, only: neko_error
8 use json_module, only: json_file
9 use json_utils, only: json_get, json_get_or_default
10 use simulation_m, only: simulation_t
11 use design, only: design_t
12 use field, only: field_t
13 use field_registry, only: neko_field_registry
14 use profiler, only: profiler_start_region, profiler_end_region
15
16 use vector, only: vector_t
17 use matrix, only: matrix_t
18
19 !only to print nglobal when running in parallel
20 use comm, only: neko_comm, pe_rank
21 use mpi_f08, only: mpi_integer, mpi_sum, mpi_allreduce
22
23 use neko_config, only: neko_bcknd_device
24 ! Inclusions from external dependencies and standard libraries
25 use, intrinsic :: iso_fortran_env, only: stderr => error_unit
26
27 use math, only: copy, cmult
28 use device_math, only: device_copy, device_cmult
29 use field_math, only: field_rzero
30 use vector_math, only: vector_cmult
31 use neko_ext, only: reset
33 use device, only: device_memcpy, host_to_device, device_to_host
34
35 implicit none
36 private
37 public :: mma_optimizer_t
38
39 ! Concrete type for MMA optimizer
40 type, extends(optimizer_t) :: mma_optimizer_t
41
42 type(mma_t) :: mma
43
50 real(kind=rp) :: scale
51 logical :: auto_scale
52
53 contains
54
55 ! Override the deferred methods
56 generic :: init => init_from_json, init_from_components
57 procedure, pass(this) :: init_from_json => mma_optimizer_init_from_json
58 procedure, pass(this) :: init_from_components => &
59 mma_optimizer_init_from_components
60
61 procedure, pass(this) :: run => mma_optimizer_run
62 procedure, pass(this) :: validate => mma_optimizer_validate
63 procedure, pass(this) :: free => mma_optimizer_free
64
65 end type mma_optimizer_t
66
67contains
68
70 subroutine mma_optimizer_init_from_json(this, parameters, problem, design, &
71 max_iterations, tolerance, simulation)
72 class(mma_optimizer_t), intent(inout) :: this
73 type(json_file), intent(inout) :: parameters
74 class(problem_t), intent(in) :: problem
75 class(design_t), intent(in) :: design
76 integer, intent(in) :: max_iterations
77 real(kind=rp), intent(in) :: tolerance
78 type(simulation_t), optional, intent(in) :: simulation
79
80 character(len=1024) :: optimization_header
81 character(len=1024) :: problem_header
82
83 type(vector_t) :: x
84 type(json_file) :: solver_parameters
85
86 ! Initialize the logger
87 call this%logger%init('optimization_data.csv')
88
89 ! Write the header
90 problem_header = problem%get_log_header()
91 optimization_header = 'iter, ' // trim(problem_header) // &
92 ', KKTmax, KKTnorm2, scaling factor'
93 call this%logger%set_header(trim(optimization_header))
94
95 call design%get_values(x)
96
97 if (pe_rank .eq. 0) then
98 print *, "Initializing mma_optimizer with steady_state_problem_t."
99 end if
100
101 call json_get(parameters, "optimization.solver", &
102 solver_parameters)
103 call this%mma%init(x%x, design%size(), problem%get_n_constraints(), &
104 solver_parameters, this%scale, this%auto_scale)
105
106 call this%init_from_components(problem, design, &
107 max_iterations, tolerance, simulation)
108
109 call x%free()
110 end subroutine mma_optimizer_init_from_json
111
113 subroutine mma_optimizer_init_from_components(this, problem, design, &
114 max_iterations, tolerance, simulation)
115 class(mma_optimizer_t), intent(inout) :: this
116 class(problem_t), intent(in) :: problem
117 class(design_t), intent(in) :: design
118 integer, intent(in) :: max_iterations
119 real(kind=rp), intent(in) :: tolerance
120 type(simulation_t), intent(in), optional :: simulation
121
122 call this%init_base(max_iterations, tolerance)
123
124 end subroutine mma_optimizer_init_from_components
125
127 subroutine mma_optimizer_run(this, problem, design, simulation)
128 class(mma_optimizer_t), intent(inout) :: this
129 class(problem_t), intent(inout) :: problem
130 class(design_t), intent(inout) :: design
131 type(simulation_t), optional, intent(inout) :: simulation
132
133 type(vector_t) :: x
134
135 integer :: iter, ierr, nglobal, n
136 real(kind=rp) :: scaling_factor
137
138 real(kind=rp) :: objective_value
139 type(vector_t) :: all_objectives
140 type(vector_t) :: constraint_value
141 type(vector_t) :: objective_sensitivities
142 type(matrix_t) :: constraint_sensitivities
143
144 type(vector_t) :: log_data
145
146 call profiler_start_region("Optimizer iteration 0")
147
148 n = design%size()
149 call mpi_allreduce(n, nglobal, 1, mpi_integer, mpi_sum, neko_comm, ierr)
150
151 ! Initialize the vectors
152 call x%init(n)
153 call all_objectives%init(problem%get_n_objectives())
154 call constraint_value%init(problem%get_n_constraints())
155 call objective_sensitivities%init(n)
156 call constraint_sensitivities%init(problem%get_n_constraints(), n)
157
159 scaling_factor = 1.0_rp
160 if (pe_rank .eq. 0) then
161 print *, "max_iterations for the optimization loop = ", &
162 this%max_iterations
163 end if
164
165 call problem%compute(design, simulation)
166 call problem%compute_sensitivity(design, simulation)
167
168 call problem%get_objective_value(objective_value)
169 call problem%get_constraint_values(constraint_value)
170 call problem%get_objective_sensitivities(objective_sensitivities)
171 call problem%get_constraint_sensitivities(constraint_sensitivities)
172 call problem%get_all_objective_values(all_objectives)
173
174 ! Stamp the initial condition
175 call mma_logger_assemble_data(log_data, 0, objective_value, &
176 all_objectives, constraint_value, 0.0_rp, 0.0_rp, scaling_factor, &
177 problem%get_n_objectives(), problem%get_n_constraints())
178 call this%logger%write(log_data)
179
180 if (present(simulation)) call simulation%write(0)
181 call design%write(0)
182 call profiler_end_region("Optimizer iteration 0")
183
184 do iter = 1, this%max_iterations
185 if (this%mma%get_residumax() .lt. this%tolerance) exit
186
187 call profiler_start_region("Optimizer iteration")
188
189 ! Scaling
190 if (this%auto_scale .eqv. .true.) then
191 scaling_factor = abs(this%scale/constraint_value%x(1))
192 else
193 scaling_factor = abs(this%scale)
194 end if
195
196 call design%get_values(x)
197
198 call vector_cmult(constraint_value, scaling_factor)
199
200 if (neko_bcknd_device .eq. 1) then
201 call device_cmult(constraint_sensitivities%x_d, scaling_factor, &
202 constraint_sensitivities%size())
203 else
204 call cmult(constraint_sensitivities%x, scaling_factor, &
205 constraint_sensitivities%size())
206 end if
207
208 ! Use scaled sensitivities to update the design variable
209 call profiler_start_region("MMA update")
210 call this%mma%update(iter, x, objective_sensitivities, &
211 constraint_value, constraint_sensitivities)
212 call profiler_end_region("MMA update")
213
214 call design%update_design(x)
215
216 call problem%compute(design, simulation)
217 call problem%compute_sensitivity(design, simulation)
218
219 call problem%get_objective_value(objective_value)
220 call problem%get_constraint_values(constraint_value)
221 call problem%get_objective_sensitivities(objective_sensitivities)
222 call problem%get_constraint_sensitivities(constraint_sensitivities)
223 call problem%get_all_objective_values(all_objectives)
224
225 call profiler_start_region("MMA KKT computation")
226 call this%mma%KKT(x, objective_sensitivities, &
227 constraint_value, constraint_sensitivities)
228 call profiler_end_region("MMA KKT computation")
229
230 ! Stamp the i^th iteration
231 call mma_logger_assemble_data(log_data, iter, objective_value, &
232 all_objectives, constraint_value, this%mma%get_residumax(), &
233 this%mma%get_residunorm(), scaling_factor, &
234 problem%get_n_objectives(), problem%get_n_constraints())
235 call this%logger%write(log_data)
236
237 if (present(simulation)) call simulation%write(iter)
238 call design%write(iter)
239
240 call profiler_end_region("Optimizer iteration")
241 end do
242
243 call this%validate(problem, design)
244
245 ! Final state after optimization
246 if (pe_rank .eq. 0) then
247 print *, "MMA Optimization completed after", iter-1, "iterations."
248 end if
249
250 ! Free local resources
251 call x%free()
252 call log_data%free()
253 call all_objectives%free()
254 call constraint_value%free()
255 call objective_sensitivities%free()
256 call constraint_sensitivities%free()
257
258 end subroutine mma_optimizer_run
259
261 subroutine mma_optimizer_validate(this, problem, design)
262 class(mma_optimizer_t), intent(inout) :: this
263 class(problem_t), intent(in) :: problem
264 class(design_t), intent(in) :: design
265
266 type(vector_t) :: constraint_values
267
268 call constraint_values%init(problem%get_n_constraints())
269 call problem%get_constraint_values(constraint_values)
270 if (neko_bcknd_device .eq. 1) then
271 call device_memcpy(constraint_values%x, constraint_values%x_d, &
272 constraint_values%size(), host_to_device, .true.)
273 end if
274
275 if (any(constraint_values%x .gt. 0.0_rp)) then
276 call neko_error("MMA optimizer validation failed: " // &
277 "Constraints are not satisfied.")
278 end if
279
280 ! Free local resources
281 call constraint_values%free()
282
283 end subroutine mma_optimizer_validate
284
285 ! Free resources associated with the MMA optimizer
286 subroutine mma_optimizer_free(this)
287 class(mma_optimizer_t), intent(inout) :: this
288
289 ! Free MMA-specific data
290 call this%mma%free()
291 end subroutine mma_optimizer_free
292
293 ! package up the log data
294 subroutine mma_logger_assemble_data(log_data, iter, objective_value, &
295 all_objectives, constraint_value, residumax, residunorm, &
296 scaling_factor, n, m)
297 type(vector_t), intent(out) :: log_data
298 integer, intent(in) :: iter
299 real(kind=rp), intent(in) ::objective_value
300 type(vector_t), intent(in) :: all_objectives
301 type(vector_t), intent(in) :: constraint_value
302 real(kind=rp), intent(in) :: residumax, residunorm, scaling_factor
303 integer, intent(in) :: n, m
304 integer :: i_tmp1, i_tmp2
305
306 ! initialize the logger data
307 ! iter | tot F | F_1 | .. |F_n | C_1 | ... | C_n | KKT | KKT2 | scale |
308 call log_data%init(5 + n + m)
309
310 ! iteration
311 log_data%x(1) = real(iter, kind=rp)
312
313 ! total objective
314 log_data%x(2) = objective_value
315
316 ! individual objectives
317 i_tmp1 = 3
318 i_tmp2 = i_tmp1 + n - 1
319 log_data%x(i_tmp1 : i_tmp2) = all_objectives%x
320
321 ! constraints
322 i_tmp1 = i_tmp2 + 1
323 i_tmp2 = i_tmp1 + m - 1
324 log_data%x(i_tmp1 : i_tmp2) = constraint_value%x
325
326 ! convergence stuff
327 log_data%x(i_tmp2 + 1) = residumax
328 log_data%x(i_tmp2 + 2) = residunorm
329 log_data%x(i_tmp2 + 3) = scaling_factor
330
331 end subroutine mma_logger_assemble_data
332end module mma_optimizer
333
Implements the design_t.
Definition design.f90:34
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 reset(neko_case)
Reset the case data structure.
Definition neko_ext.f90:58
Module for handling the optimization problem.
Definition problem.f90:35
Implements the steady_problem_t type.
An abstract design type.
Definition design.f90:52
Abstract optimizer class.
Definition optimizer.f90:18
The abstract problem type.
Definition problem.f90:61