42 use num_types,
only: rp, dp
47 use vector,
only: vector_t
48 use matrix,
only: matrix_t
49 use device,
only: host_to_device, device_to_host
50 use json_module,
only: json_file
51 use json_utils,
only: json_extract_item, json_get, json_get_or_default
53 use logger,
only: neko_log
55 use vector_math,
only: vector_add2, vector_cfill
56 use time_step_controller,
only: time_step_controller_t
59 use simulation,
only: simulation_init, simulation_step, simulation_finalize
60 use mpi_f08,
only: mpi_wtime
61 use profiler,
only: profiler_start_region, profiler_end_region
70 integer :: n_design = 0
72 integer :: n_objectives = 0
74 integer :: n_constraints = 0
88 procedure, pass(this),
public :: free => problem_free
93 procedure, pass(this),
public :: compute => problem_compute
98 procedure, pass(this),
public :: compute_sensitivity => &
99 problem_compute_sensitivity
102 procedure, pass(this),
public :: run_forward_unsteady => &
103 problem_run_forward_unsteady
106 procedure, pass(this),
public :: run_backward_unsteady => &
107 problem_run_backward_unsteady
112 procedure, pass(this),
public :: read_objectives => problem_read_objectives
114 procedure, pass(this),
public :: read_constraints => &
115 problem_read_constraints
121 procedure, pass(this),
public :: write => problem_write
123 procedure, pass(this),
public :: get_log_values => problem_get_log_values
126 procedure, pass(this),
public :: add_objective => problem_add_objective
128 procedure, pass(this),
public :: add_constraint => problem_add_constraint
134 procedure, pass(this) :: update_objectives => &
135 problem_update_objectives
137 procedure, pass(this) :: update_constraints => &
138 problem_update_constraints
140 procedure, pass(this) :: update_objective_sensitivities => &
141 problem_update_objective_sensitivities
143 procedure, pass(this) :: update_constraint_sensitivities => &
144 problem_update_constraint_sensitivities
147 procedure, pass(this) :: reset_objectives => &
148 problem_reset_objectives
150 procedure, pass(this) :: reset_constraints => &
151 problem_reset_constraints
153 procedure, pass(this) :: reset_objective_sensitivities => &
154 problem_reset_objective_sensitivities
156 procedure, pass(this) :: reset_constraint_sensitivities => &
157 problem_reset_constraint_sensitivities
160 procedure, pass(this) :: accumulate_objectives => &
161 problem_accumulate_objectives
163 procedure, pass(this) :: accumulate_constraints => &
164 problem_accumulate_constraints
166 procedure, pass(this) :: accumulate_objective_sensitivities => &
167 problem_accumulate_objective_sensitivities
169 procedure, pass(this) :: accumulate_constraint_sensitivities => &
170 problem_accumulate_constraint_sensitivities
176 procedure, pass(this),
public :: get_objective_value => &
177 problem_get_objective_value
179 procedure, pass(this),
public :: get_all_objective_values => &
180 problem_get_all_objective_values
182 procedure, pass(this),
public :: get_constraint_values => &
183 problem_get_constraint_values
185 procedure, pass(this),
public :: get_objective_sensitivities => &
186 problem_get_objective_sensitivities
188 procedure, pass(this),
public :: get_constraint_sensitivities => &
189 problem_get_constraint_sensitivities
192 procedure, pass(this) :: get_n_objectives => problem_get_num_objectives
194 procedure, pass(this) :: get_n_constraints => problem_get_num_constraints
197 procedure, pass(this) :: get_log_header => problem_get_log_header
199 procedure, pass(this) :: get_log_size => problem_get_log_size
211 type(json_file),
intent(inout) :: parameters
212 class(
design_t),
intent(in) :: design
213 type(
simulation_t),
optional,
intent(inout) :: simulation
217 this%n_design =
design%size()
220 call this%read_objectives(parameters,
design, simulation)
221 call this%read_constraints(parameters,
design, simulation)
226 subroutine problem_free(this)
231 this%n_objectives = 0
232 this%n_constraints = 0
235 if (
allocated(this%objective_list))
then
236 do i = 1,
size(this%objective_list)
237 call this%objective_list(i)%free()
239 deallocate(this%objective_list)
243 if (
allocated(this%constraint_list))
then
244 do i = 1,
size(this%constraint_list)
245 call this%constraint_list(i)%free()
247 deallocate(this%constraint_list)
249 end subroutine problem_free
252 subroutine problem_write(this, idx)
254 integer,
intent(in) :: idx
256 end subroutine problem_write
262 subroutine problem_get_log_values(this, values, include_constraints)
264 real(kind=rp),
intent(out) :: values(:)
265 logical,
intent(in),
optional :: include_constraints
266 integer :: i, n, offset
267 real(kind=rp) :: objective_value
268 real(kind=rp),
allocatable :: tmp(:)
269 logical :: do_constraints
271 if (
present(include_constraints))
then
272 do_constraints = include_constraints
274 do_constraints = .true.
277 call this%get_objective_value(objective_value)
279 values(1) = objective_value
282 do i = 1, this%n_objectives
283 n = this%objective_list(i)%objective%get_log_size()
286 call this%objective_list(i)%objective%get_log_values(tmp)
287 values(offset:offset + n - 1) = tmp
293 if (do_constraints)
then
294 do i = 1, this%n_constraints
295 n = this%constraint_list(i)%constraint%get_log_size()
298 call this%constraint_list(i)%constraint%get_log_values(tmp)
299 values(offset:offset + n - 1) = tmp
305 end subroutine problem_get_log_values
311 subroutine problem_read_objectives(this, parameters, design, simulation)
313 type(json_file),
intent(inout) :: parameters
314 class(
design_t),
intent(in) :: design
315 type(
simulation_t),
optional,
intent(inout) :: simulation
319 character(len=:),
allocatable :: path, type
320 type(json_file) :: objective_json
321 integer :: n_objectives, i
324 call neko_log%section(
"Reading objectives")
327 path =
"optimization.objectives"
328 if (parameters%valid_path(path))
then
329 call parameters%info(path, n_children = n_objectives)
332 do i = 1, n_objectives
333 call json_extract_item(parameters, path, i, objective_json)
334 call json_get(objective_json,
"type", type)
335 call neko_log%message(type)
342 if (
present(simulation))
then
347 call json_get_or_default(parameters, &
348 "adjoint_fluid.dealias_sensitivity", dealias, .true.)
349 call alo%init_from_attributes(
design, simulation, weight = 1.0_rp, &
350 name =
"Augmented Lagrangian", mask_name =
"", &
356 call neko_log%end_section()
358 end subroutine problem_read_objectives
361 subroutine problem_read_constraints(this, parameters, design, simulation)
363 type(json_file),
intent(inout) :: parameters
364 class(
design_t),
intent(in) :: design
366 type(
simulation_t),
optional,
intent(inout) :: simulation
369 character(len=:),
allocatable :: path, type
370 type(json_file) :: constraint_json
371 integer :: n_constraints, i
373 call neko_log%section(
"Reading constraints")
376 path =
"optimization.constraints"
378 if (parameters%valid_path(path))
then
379 call parameters%info(path, n_children = n_constraints)
382 do i = 1, n_constraints
383 call json_extract_item(parameters, path, i, constraint_json)
384 call json_get(constraint_json,
"type", type)
385 call neko_log%message(type)
392 call neko_log%end_section()
394 end subroutine problem_read_constraints
397 subroutine problem_add_objective(this, objective)
399 class(
objective_t),
allocatable,
intent(inout) :: objective
404 if (
allocated(this%objective_list))
then
405 n =
size(this%objective_list)
406 call move_alloc(this%objective_list, temp_list)
407 allocate(this%objective_list(n + 1))
408 if (
allocated(temp_list))
then
410 call move_alloc(temp_list(i)%objective, &
411 this%objective_list(i)%objective)
415 allocate(this%objective_list(1))
418 call move_alloc(
objective, this%objective_list(n + 1)%objective)
419 this%n_objectives = n + 1
420 end subroutine problem_add_objective
423 subroutine problem_add_constraint(this, constraint)
425 class(
constraint_t),
allocatable,
intent(inout) :: constraint
430 if (
allocated(this%constraint_list))
then
431 n =
size(this%constraint_list)
432 call move_alloc(this%constraint_list, temp_list)
433 allocate(this%constraint_list(n + 1))
434 if (
allocated(temp_list))
then
436 call move_alloc(temp_list(i)%constraint, &
437 this%constraint_list(i)%constraint)
441 allocate(this%constraint_list(1))
444 call move_alloc(
constraint, this%constraint_list(n + 1)%constraint)
445 this%n_constraints = n + 1
446 end subroutine problem_add_constraint
452 subroutine problem_compute(this, design, simulation)
454 class(
design_t),
intent(inout) :: design
455 class(
simulation_t),
optional,
intent(inout) :: simulation
457 if (
present(simulation))
then
458 call simulation%reset()
459 if (simulation%unsteady)
then
461 call this%run_forward_unsteady(simulation,
design)
463 call simulation%run_forward()
465 call this%update_objectives(
design)
468 call this%update_objectives(
design)
471 call this%update_constraints(
design)
473 end subroutine problem_compute
476 subroutine problem_compute_sensitivity(this, design, simulation)
478 class(
design_t),
intent(inout) :: design
479 class(
simulation_t),
optional,
intent(inout) :: simulation
481 type(vector_t) :: objective_sensitivity
483 if (
present(simulation))
then
484 if (simulation%unsteady)
then
486 call this%run_backward_unsteady(simulation,
design)
488 call simulation%run_backward()
490 call this%update_objective_sensitivities(
design)
493 call this%update_objective_sensitivities(
design)
496 call this%update_constraint_sensitivities(
design)
498 call objective_sensitivity%init(this%n_design)
499 call this%get_objective_sensitivities(objective_sensitivity)
501 call design%map_backward(objective_sensitivity)
503 call objective_sensitivity%free()
504 end subroutine problem_compute_sensitivity
507 subroutine problem_run_forward_unsteady(this, simulation, design)
510 class(
design_t),
intent(inout) :: design
511 type(time_step_controller_t) :: dt_controller
512 real(kind=dp) :: loop_start
514 call dt_controller%init(simulation%neko_case%params)
516 call simulation%reset()
517 call simulation_init(simulation%neko_case, dt_controller)
520 call this%reset_objectives()
522 call profiler_start_region(
"Forward simulation")
523 loop_start = mpi_wtime()
524 simulation%n_timesteps = 0
525 do while (simulation%neko_case%time%t .lt. simulation%neko_case%time%end_time)
526 simulation%n_timesteps = simulation%n_timesteps + 1
528 call simulation_step(simulation%neko_case, dt_controller, loop_start)
530 call this%accumulate_objectives(
design, simulation%adjoint_case%time%dt)
532 call simulation%checkpoint%save(simulation%neko_case)
534 call profiler_end_region(
"Forward simulation")
536 call simulation_finalize(simulation%neko_case)
538 end subroutine problem_run_forward_unsteady
541 subroutine problem_run_backward_unsteady(this, simulation, design)
544 class(
design_t),
intent(inout) :: design
545 type(time_step_controller_t) :: dt_controller
546 real(kind=dp) :: loop_start
548 real(kind=rp) :: total_time
551 call dt_controller%init(simulation%neko_case%params)
556 call this%reset_objective_sensitivities()
558 cfl = simulation%adjoint_case%fluid_adj%compute_cfl(simulation%adjoint_case%time%dt)
559 loop_start = mpi_wtime()
561 total_time = simulation%n_timesteps * simulation%adjoint_case%time%dt
563 call profiler_start_region(
"Adjoint simulation")
582 do i = simulation%n_timesteps, 1, -1
584 call simulation%checkpoint%restore(simulation%neko_case, i)
586 call this%accumulate_objective_sensitivities(
design, &
587 simulation%adjoint_case%time%dt)
590 cfl, loop_start, total_time)
593 call profiler_end_region(
"Forward simulation")
597 end subroutine problem_run_backward_unsteady
608 subroutine problem_update_objectives(this, design)
610 class(
design_t),
intent(in) :: design
613 do i = 1, this%n_objectives
614 call this%objective_list(i)%objective%update_value(
design)
616 end subroutine problem_update_objectives
624 subroutine problem_update_constraints(this, design)
626 class(
design_t),
intent(in) :: design
629 do i = 1, this%n_constraints
630 call this%constraint_list(i)%constraint%update_value(
design)
632 end subroutine problem_update_constraints
640 subroutine problem_update_objective_sensitivities(this, design)
642 class(
design_t),
intent(in) :: design
645 do i = 1, this%n_objectives
646 call this%objective_list(i)%objective%update_sensitivity(
design)
648 end subroutine problem_update_objective_sensitivities
656 subroutine problem_update_constraint_sensitivities(this, design)
658 class(
design_t),
intent(in) :: design
661 do i = 1, this%n_constraints
662 call this%constraint_list(i)%constraint%update_sensitivity(
design)
664 end subroutine problem_update_constraint_sensitivities
673 subroutine problem_reset_objectives(this)
677 do i = 1, this%n_objectives
678 call this%objective_list(i)%objective%reset_value()
680 end subroutine problem_reset_objectives
686 subroutine problem_reset_constraints(this)
690 do i = 1, this%n_constraints
691 call this%constraint_list(i)%constraint%reset_value()
693 end subroutine problem_reset_constraints
699 subroutine problem_reset_objective_sensitivities(this)
703 do i = 1, this%n_objectives
704 call this%objective_list(i)%objective%reset_sensitivity()
706 end subroutine problem_reset_objective_sensitivities
712 subroutine problem_reset_constraint_sensitivities(this)
716 do i = 1, this%n_constraints
717 call this%constraint_list(i)%constraint%reset_sensitivity()
719 end subroutine problem_reset_constraint_sensitivities
730 subroutine problem_accumulate_objectives(this, design, dt)
732 class(
design_t),
intent(in) :: design
733 real(kind=rp),
intent(in) :: dt
736 do i = 1, this%n_objectives
737 call this%objective_list(i)%objective%accumulate_value(
design, dt)
739 end subroutine problem_accumulate_objectives
747 subroutine problem_accumulate_constraints(this, design, dt)
749 class(
design_t),
intent(in) :: design
750 real(kind=rp),
intent(in) :: dt
753 do i = 1, this%n_constraints
754 call this%constraint_list(i)%constraint%accumulate_value(
design, dt)
756 end subroutine problem_accumulate_constraints
764 subroutine problem_accumulate_objective_sensitivities(this, design, dt)
766 class(
design_t),
intent(in) :: design
767 real(kind=rp),
intent(in) :: dt
770 do i = 1, this%n_objectives
771 call this%objective_list(i)%objective%accumulate_sensitivity(
design, dt)
773 end subroutine problem_accumulate_objective_sensitivities
781 subroutine problem_accumulate_constraint_sensitivities(this, design, dt)
783 class(
design_t),
intent(in) :: design
784 real(kind=rp),
intent(in) :: dt
787 do i = 1, this%n_constraints
788 call this%constraint_list(i)%constraint%accumulate_sensitivity(
design, dt)
790 end subroutine problem_accumulate_constraint_sensitivities
801 subroutine problem_get_objective_value(this, objective_value)
803 real(kind=rp),
intent(out) :: objective_value
806 objective_value = 0.0_rp
807 do i = 1, this%n_objectives
808 objective_value = objective_value + &
809 this%objective_list(i)%objective%get_weight() * &
810 this%objective_list(i)%objective%get_value()
813 end subroutine problem_get_objective_value
821 subroutine problem_get_all_objective_values(this, all_objective_values)
823 type(vector_t),
intent(inout) :: all_objective_values
826 do i = 1, this%n_objectives
827 all_objective_values%x(i) = this%objective_list(i)%objective%value
830 call all_objective_values%copy_from(host_to_device, sync = .true.)
832 end subroutine problem_get_all_objective_values
840 subroutine problem_get_constraint_values(this, constraint_value)
842 type(vector_t),
intent(inout) :: constraint_value
845 do i = 1, this%n_constraints
846 constraint_value%x(i) = this%constraint_list(i)%constraint%value
849 call constraint_value%copy_from(host_to_device, sync = .true.)
851 end subroutine problem_get_constraint_values
859 subroutine problem_get_objective_sensitivities(this, sensitivity)
861 type(vector_t),
intent(inout) :: sensitivity
864 call vector_cfill(sensitivity, 0.0_rp)
865 do i = 1, this%n_objectives
866 call vector_add2(sensitivity, &
867 this%objective_list(i)%objective%sensitivity)
870 end subroutine problem_get_objective_sensitivities
878 subroutine problem_get_constraint_sensitivities(this, sensitivity)
880 type(matrix_t),
target,
intent(inout) :: sensitivity
881 real(kind=rp),
pointer :: row(:)
885 do i = 1, this%n_constraints
886 call this%constraint_list(i)%constraint%sensitivity%copy_from( &
887 device_to_host, sync = i .eq. this%n_constraints)
890 do i = 1, this%n_constraints
891 row(1:this%n_design) => sensitivity%x(i, :)
893 call copy(row, this%constraint_list(i)%constraint%sensitivity%x, &
897 call sensitivity%copy_from(host_to_device, sync = .true.)
899 end subroutine problem_get_constraint_sensitivities
905 pure function problem_get_num_objectives(this)
result(n)
909 n = this%n_objectives
910 end function problem_get_num_objectives
913 pure function problem_get_num_constraints(this)
result(n)
917 n = this%n_constraints
918 end function problem_get_num_constraints
924 function problem_get_log_header(this, include_constraints)
result(buff)
926 logical,
intent(in),
optional :: include_constraints
927 character(len=4096) :: buff
928 character(len=128) :: mini_buff
929 character(len=128),
allocatable :: headers(:)
931 logical :: do_constraints
933 buff =
"Total objective function"
934 if (
present(include_constraints))
then
935 do_constraints = include_constraints
937 do_constraints = .true.
940 do i = 1, this%get_n_objectives()
941 n = this%objective_list(i)%objective%get_log_size()
944 call this%objective_list(i)%objective%get_log_headers(headers)
947 write(mini_buff,
'(", ", A)') trim(headers(j))
948 buff = trim(buff) // trim(mini_buff)
954 if (do_constraints)
then
955 do i = 1, this%get_n_constraints()
956 n = this%constraint_list(i)%constraint%get_log_size()
959 call this%constraint_list(i)%constraint%get_log_headers(headers)
962 write(mini_buff,
'(", ", A)') trim(headers(j))
963 buff = trim(buff) // trim(mini_buff)
970 end function problem_get_log_header
976 function problem_get_log_size(this, include_constraints)
result(n)
978 logical,
intent(in),
optional :: include_constraints
980 logical :: do_constraints
983 if (
present(include_constraints))
then
984 do_constraints = include_constraints
986 do_constraints = .true.
988 do i = 1, this%get_n_objectives()
989 n = n + this%objective_list(i)%objective%get_log_size()
992 if (do_constraints)
then
993 do i = 1, this%get_n_constraints()
994 n = n + this%constraint_list(i)%constraint%get_log_size()
998 end function problem_get_log_size
Factory function Allocates and initializes an constraint function object.
Factory function Allocates and initializes an objective function object.
Implements the augmented_lagrangian_objective_t type.
Implements the constraint_t type.
Implements the objective_t type.
Module for handling the optimization problem.
subroutine problem_init(this, parameters, design, simulation)
The constructor for the base problem.
Adjoint simulation driver.
subroutine, public simulation_adjoint_init(c, dt_controller)
Initialise a simulation_adjoint of a case.
subroutine, public simulation_adjoint_step(c, dt_controller, cfl, tstep_loop_start_time, final_time)
Compute a single time-step of an adjoint case.
subroutine, public simulation_adjoint_finalize(c)
Finalize a simulation of a case.
Implements the steady_problem_t type.
An objective function implementing our augmented lagrangian sensitivity contribution.
The abstract constraint type.
Wrapper for constraints for use in lists.
The abstract objective type.
Wrapper for objectives for use in lists.
The abstract problem type.