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
12 use field,
only: field_t
13 use field_registry,
only: neko_field_registry
14 use profiler,
only: profiler_start_region, profiler_end_region
16 use vector,
only: vector_t
17 use matrix,
only: matrix_t
20 use comm,
only: neko_comm, pe_rank
21 use mpi_f08,
only: mpi_integer, mpi_sum, mpi_allreduce
23 use neko_config,
only: neko_bcknd_device
25 use,
intrinsic :: iso_fortran_env, only: stderr => error_unit
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
33 use device,
only: device_memcpy, host_to_device, device_to_host
50 real(kind=rp) :: scale
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
61 procedure, pass(this) :: run => mma_optimizer_run
62 procedure, pass(this) :: validate => mma_optimizer_validate
63 procedure, pass(this) :: free => mma_optimizer_free
70 subroutine mma_optimizer_init_from_json(this, parameters, problem, design, &
71 max_iterations, tolerance, simulation)
73 type(json_file),
intent(inout) :: parameters
76 integer,
intent(in) :: max_iterations
77 real(kind=rp),
intent(in) :: tolerance
80 character(len=1024) :: optimization_header
81 character(len=1024) :: problem_header
84 type(json_file) :: solver_parameters
87 call this%logger%init(
'optimization_data.csv')
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))
97 if (pe_rank .eq. 0)
then
98 print *,
"Initializing mma_optimizer with steady_state_problem_t."
101 call json_get(parameters,
"optimization.solver", &
103 call this%mma%init(x%x,
design%size(),
problem%get_n_constraints(), &
104 solver_parameters, this%scale, this%auto_scale)
107 max_iterations, tolerance, simulation)
110 end subroutine mma_optimizer_init_from_json
113 subroutine mma_optimizer_init_from_components(this, problem, design, &
114 max_iterations, tolerance, simulation)
118 integer,
intent(in) :: max_iterations
119 real(kind=rp),
intent(in) :: tolerance
122 call this%init_base(max_iterations, tolerance)
124 end subroutine mma_optimizer_init_from_components
127 subroutine mma_optimizer_run(this, problem, design, simulation)
131 type(
simulation_t),
optional,
intent(inout) :: simulation
135 integer :: iter, ierr, nglobal, n
136 real(kind=rp) :: scaling_factor
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
144 type(vector_t) :: log_data
146 call profiler_start_region(
"Optimizer iteration 0")
149 call mpi_allreduce(n, nglobal, 1, mpi_integer, mpi_sum, neko_comm, ierr)
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)
159 scaling_factor = 1.0_rp
160 if (pe_rank .eq. 0)
then
161 print *,
"max_iterations for the optimization loop = ", &
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)
175 call mma_logger_assemble_data(log_data, 0, objective_value, &
176 all_objectives, constraint_value, 0.0_rp, 0.0_rp, scaling_factor, &
178 call this%logger%write(log_data)
180 if (
present(simulation))
call simulation%write(0)
182 call profiler_end_region(
"Optimizer iteration 0")
184 do iter = 1, this%max_iterations
185 if (this%mma%get_residumax() .lt. this%tolerance)
exit
187 call profiler_start_region(
"Optimizer iteration")
190 if (this%auto_scale .eqv. .true.)
then
191 scaling_factor = abs(this%scale/constraint_value%x(1))
193 scaling_factor = abs(this%scale)
198 call vector_cmult(constraint_value, scaling_factor)
200 if (neko_bcknd_device .eq. 1)
then
201 call device_cmult(constraint_sensitivities%x_d, scaling_factor, &
202 constraint_sensitivities%size())
204 call cmult(constraint_sensitivities%x, scaling_factor, &
205 constraint_sensitivities%size())
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")
214 call design%update_design(x)
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)
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")
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, &
235 call this%logger%write(log_data)
237 if (
present(simulation))
call simulation%write(iter)
240 call profiler_end_region(
"Optimizer iteration")
246 if (pe_rank .eq. 0)
then
247 print *,
"MMA Optimization completed after", iter-1,
"iterations."
253 call all_objectives%free()
254 call constraint_value%free()
255 call objective_sensitivities%free()
256 call constraint_sensitivities%free()
258 end subroutine mma_optimizer_run
261 subroutine mma_optimizer_validate(this, problem, design)
266 type(vector_t) :: constraint_values
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.)
275 if (any(constraint_values%x .gt. 0.0_rp))
then
276 call neko_error(
"MMA optimizer validation failed: " // &
277 "Constraints are not satisfied.")
281 call constraint_values%free()
283 end subroutine mma_optimizer_validate
286 subroutine mma_optimizer_free(this)
291 end subroutine mma_optimizer_free
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
308 call log_data%init(5 + n + m)
311 log_data%x(1) = real(iter, kind=rp)
314 log_data%x(2) = objective_value
318 i_tmp2 = i_tmp1 + n - 1
319 log_data%x(i_tmp1 : i_tmp2) = all_objectives%x
323 i_tmp2 = i_tmp1 + m - 1
324 log_data%x(i_tmp1 : i_tmp2) = constraint_value%x
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
331 end subroutine mma_logger_assemble_data
332end module mma_optimizer
Some common Masking operations we may need.
Contains extensions to the neko library required to run the topology optimization code.
subroutine, public reset(neko_case)
Reset the case data structure.
Module for handling the optimization problem.
Implements the steady_problem_t type.
Abstract optimizer class.
The abstract problem type.