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 time_state,
only: time_state_t
56 use vector_math,
only: vector_add2, vector_cfill
57 use time_step_controller,
only: time_step_controller_t
60 use simulation,
only: simulation_init, simulation_step, simulation_finalize
61 use mpi_f08,
only: mpi_wtime
62 use profiler,
only: profiler_start_region, profiler_end_region
63 use utils,
only: neko_error
72 integer :: n_design = 0
74 integer :: n_objectives = 0
76 integer :: n_constraints = 0
90 procedure, pass(this),
public :: free => problem_free
95 procedure, pass(this),
public :: compute => problem_compute
100 procedure, pass(this),
public :: compute_sensitivity => &
101 problem_compute_sensitivity
104 procedure, pass(this),
public :: run_forward_unsteady => &
105 problem_run_forward_unsteady
108 procedure, pass(this),
public :: run_backward_unsteady => &
109 problem_run_backward_unsteady
114 procedure, pass(this),
public :: read_objectives => problem_read_objectives
116 procedure, pass(this),
public :: read_constraints => &
117 problem_read_constraints
123 procedure, pass(this),
public :: write => problem_write
125 procedure, pass(this),
public :: get_log_values => problem_get_log_values
128 procedure, pass(this),
public :: add_objective => problem_add_objective
130 procedure, pass(this),
public :: add_constraint => problem_add_constraint
136 procedure, pass(this) :: update_objectives => &
137 problem_update_objectives
139 procedure, pass(this) :: update_constraints => &
140 problem_update_constraints
142 procedure, pass(this) :: update_objective_sensitivities => &
143 problem_update_objective_sensitivities
145 procedure, pass(this) :: update_constraint_sensitivities => &
146 problem_update_constraint_sensitivities
149 procedure, pass(this) :: reset_objectives => &
150 problem_reset_objectives
152 procedure, pass(this) :: reset_constraints => &
153 problem_reset_constraints
155 procedure, pass(this) :: reset_objective_sensitivities => &
156 problem_reset_objective_sensitivities
158 procedure, pass(this) :: reset_constraint_sensitivities => &
159 problem_reset_constraint_sensitivities
162 procedure, pass(this) :: accumulate_objectives => &
163 problem_accumulate_objectives
165 procedure, pass(this) :: accumulate_constraints => &
166 problem_accumulate_constraints
168 procedure, pass(this) :: accumulate_objective_sensitivities => &
169 problem_accumulate_objective_sensitivities
171 procedure, pass(this) :: accumulate_constraint_sensitivities => &
172 problem_accumulate_constraint_sensitivities
178 procedure, pass(this),
public :: get_objective_value => &
179 problem_get_objective_value
181 procedure, pass(this),
public :: get_all_objective_values => &
182 problem_get_all_objective_values
184 procedure, pass(this),
public :: get_constraint_values => &
185 problem_get_constraint_values
187 procedure, pass(this),
public :: get_objective_sensitivities => &
188 problem_get_objective_sensitivities
190 procedure, pass(this),
public :: get_constraint_sensitivities => &
191 problem_get_constraint_sensitivities
194 procedure, pass(this) :: get_n_objectives => problem_get_num_objectives
196 procedure, pass(this) :: get_n_constraints => problem_get_num_constraints
199 procedure, pass(this) :: get_log_header => problem_get_log_header
201 procedure, pass(this) :: get_log_size => problem_get_log_size
213 type(json_file),
intent(inout) :: parameters
214 class(
design_t),
intent(in) :: design
215 type(
simulation_t),
optional,
intent(inout) :: simulation
219 this%n_design =
design%size()
222 call this%read_objectives(parameters,
design, simulation)
223 call this%read_constraints(parameters,
design, simulation)
228 subroutine problem_free(this)
233 this%n_objectives = 0
234 this%n_constraints = 0
237 if (
allocated(this%objective_list))
then
238 do i = 1,
size(this%objective_list)
239 call this%objective_list(i)%free()
241 deallocate(this%objective_list)
245 if (
allocated(this%constraint_list))
then
246 do i = 1,
size(this%constraint_list)
247 call this%constraint_list(i)%free()
249 deallocate(this%constraint_list)
251 end subroutine problem_free
254 subroutine problem_write(this, idx)
256 integer,
intent(in) :: idx
258 end subroutine problem_write
264 subroutine problem_get_log_values(this, values, include_constraints)
266 real(kind=rp),
intent(out) :: values(:)
267 logical,
intent(in),
optional :: include_constraints
268 integer :: i, n, offset
269 real(kind=rp) :: objective_value
270 real(kind=rp),
allocatable :: tmp(:)
271 logical :: do_constraints
273 if (
present(include_constraints))
then
274 do_constraints = include_constraints
276 do_constraints = .true.
279 call this%get_objective_value(objective_value)
281 values(1) = objective_value
284 do i = 1, this%n_objectives
285 n = this%objective_list(i)%objective%get_log_size()
288 call this%objective_list(i)%objective%get_log_values(tmp)
289 values(offset:offset + n - 1) = tmp
295 if (do_constraints)
then
296 do i = 1, this%n_constraints
297 n = this%constraint_list(i)%constraint%get_log_size()
300 call this%constraint_list(i)%constraint%get_log_values(tmp)
301 values(offset:offset + n - 1) = tmp
307 end subroutine problem_get_log_values
313 subroutine problem_read_objectives(this, parameters, design, simulation)
315 type(json_file),
intent(inout) :: parameters
316 class(
design_t),
intent(in) :: design
317 type(
simulation_t),
optional,
intent(inout) :: simulation
321 character(len=:),
allocatable :: path, type
322 type(json_file) :: objective_json
323 integer :: n_objectives, i
326 call neko_log%section(
"Reading objectives")
329 path =
"optimization.objectives"
330 if (parameters%valid_path(path))
then
331 call parameters%info(path, n_children = n_objectives)
334 do i = 1, n_objectives
335 call json_extract_item(parameters, path, i, objective_json)
336 call json_get(objective_json,
"type", type)
337 call neko_log%message(type)
344 if (
present(simulation))
then
349 call json_get_or_default(parameters, &
350 "adjoint_fluid.dealias_sensitivity", dealias, .true.)
351 call alo%init_from_attributes(
design, simulation, weight = 1.0_rp, &
352 name =
"Augmented Lagrangian", mask_name =
"", &
358 call neko_log%end_section()
360 end subroutine problem_read_objectives
363 subroutine problem_read_constraints(this, parameters, design, simulation)
365 type(json_file),
intent(inout) :: parameters
366 class(
design_t),
intent(in) :: design
368 type(
simulation_t),
optional,
intent(inout) :: simulation
371 character(len=:),
allocatable :: path, type
372 type(json_file) :: constraint_json
373 integer :: n_constraints, i
375 call neko_log%section(
"Reading constraints")
378 path =
"optimization.constraints"
380 if (parameters%valid_path(path))
then
381 call parameters%info(path, n_children = n_constraints)
384 do i = 1, n_constraints
385 call json_extract_item(parameters, path, i, constraint_json)
386 call json_get(constraint_json,
"type", type)
387 call neko_log%message(type)
395 call neko_log%end_section()
397 end subroutine problem_read_constraints
400 subroutine problem_add_objective(this, objective)
402 class(
objective_t),
allocatable,
intent(inout) :: objective
407 if (
allocated(this%objective_list))
then
408 n =
size(this%objective_list)
409 call move_alloc(this%objective_list, temp_list)
410 allocate(this%objective_list(n + 1))
411 if (
allocated(temp_list))
then
413 call move_alloc(temp_list(i)%objective, &
414 this%objective_list(i)%objective)
418 allocate(this%objective_list(1))
421 call move_alloc(
objective, this%objective_list(n + 1)%objective)
422 this%n_objectives = n + 1
423 end subroutine problem_add_objective
426 subroutine problem_add_constraint(this, constraint)
428 class(
constraint_t),
allocatable,
intent(inout) :: constraint
433 if (
allocated(this%constraint_list))
then
434 n =
size(this%constraint_list)
435 call move_alloc(this%constraint_list, temp_list)
436 allocate(this%constraint_list(n + 1))
437 if (
allocated(temp_list))
then
439 call move_alloc(temp_list(i)%constraint, &
440 this%constraint_list(i)%constraint)
444 allocate(this%constraint_list(1))
447 call move_alloc(
constraint, this%constraint_list(n + 1)%constraint)
448 this%n_constraints = n + 1
449 end subroutine problem_add_constraint
455 subroutine problem_compute(this, design, simulation)
457 class(
design_t),
intent(inout) :: design
458 class(
simulation_t),
optional,
intent(inout) :: simulation
460 if (
present(simulation))
then
461 call simulation%reset()
462 if (simulation%unsteady)
then
464 call this%run_forward_unsteady(simulation,
design)
466 call simulation%run_forward()
468 call this%update_objectives(
design)
471 call this%update_objectives(
design)
474 call this%update_constraints(
design)
476 end subroutine problem_compute
479 subroutine problem_compute_sensitivity(this, design, simulation)
481 class(
design_t),
intent(inout) :: design
482 class(
simulation_t),
optional,
intent(inout) :: simulation
484 type(vector_t) :: objective_sensitivity
486 if (
present(simulation))
then
487 if (simulation%unsteady)
then
489 call this%run_backward_unsteady(simulation,
design)
491 call simulation%run_backward()
493 call this%update_objective_sensitivities(
design)
496 call this%update_objective_sensitivities(
design)
499 call this%update_constraint_sensitivities(
design)
501 call objective_sensitivity%init(this%n_design)
502 call this%get_objective_sensitivities(objective_sensitivity)
504 call design%map_backward(objective_sensitivity)
506 call objective_sensitivity%free()
507 end subroutine problem_compute_sensitivity
510 subroutine problem_run_forward_unsteady(this, simulation, design)
513 class(
design_t),
intent(inout) :: design
514 type(time_step_controller_t) :: dt_controller
515 real(kind=dp) :: loop_start
517 call dt_controller%init(simulation%neko_case%params)
519 call simulation%reset()
520 call simulation_init(simulation%neko_case, dt_controller)
523 call this%reset_objectives()
525 if (.not.
allocated(simulation%state_recover))
then
526 call neko_error(
"State recovery not initialized.")
529 call profiler_start_region(
"Forward simulation")
530 loop_start = mpi_wtime()
531 simulation%n_timesteps = 0
532 do while (simulation%neko_case%time%t .lt. &
533 simulation%neko_case%time%end_time)
534 simulation%n_timesteps = simulation%n_timesteps + 1
536 call simulation_step(simulation%neko_case, dt_controller, loop_start)
538 call this%accumulate_objectives(
design, simulation%neko_case%time)
540 call simulation%state_recover%save()
542 call profiler_end_region(
"Forward simulation")
544 call simulation_finalize(simulation%neko_case)
546 end subroutine problem_run_forward_unsteady
549 subroutine problem_run_backward_unsteady(this, simulation, design)
552 class(
design_t),
intent(inout) :: design
553 type(time_step_controller_t) :: dt_controller
554 real(kind=dp) :: loop_start
556 real(kind=rp) :: total_time
558 type(time_state_t) :: accumulation_time
560 call dt_controller%init(simulation%neko_case%params)
565 call this%reset_objective_sensitivities()
567 cfl = simulation%adjoint_case%fluid_adj%compute_cfl( &
568 simulation%adjoint_case%time%dt)
569 loop_start = mpi_wtime()
571 if (.not.
allocated(simulation%state_recover))
then
572 call neko_error(
"State recovery not initialized.")
576 total_time = simulation%n_timesteps * simulation%adjoint_case%time%dt
578 call profiler_start_region(
"Adjoint simulation")
580 do i = simulation%n_timesteps, 1, -1
582 call simulation%state_recover%restore(i)
584 accumulation_time = simulation%adjoint_case%time
585 accumulation_time%t = total_time - simulation%adjoint_case%time%t
586 call this%accumulate_objective_sensitivities(
design, accumulation_time)
589 cfl, loop_start, total_time)
592 call profiler_end_region(
"Adjoint simulation")
596 end subroutine problem_run_backward_unsteady
607 subroutine problem_update_objectives(this, design)
609 class(
design_t),
intent(in) :: design
612 do i = 1, this%n_objectives
613 call this%objective_list(i)%objective%update_value(
design)
615 end subroutine problem_update_objectives
623 subroutine problem_update_constraints(this, design)
625 class(
design_t),
intent(in) :: design
628 do i = 1, this%n_constraints
629 call this%constraint_list(i)%constraint%update_value(
design)
631 end subroutine problem_update_constraints
639 subroutine problem_update_objective_sensitivities(this, design)
641 class(
design_t),
intent(in) :: design
644 do i = 1, this%n_objectives
645 call this%objective_list(i)%objective%update_sensitivity(
design)
647 end subroutine problem_update_objective_sensitivities
655 subroutine problem_update_constraint_sensitivities(this, design)
657 class(
design_t),
intent(in) :: design
660 do i = 1, this%n_constraints
661 call this%constraint_list(i)%constraint%update_sensitivity(
design)
663 end subroutine problem_update_constraint_sensitivities
672 subroutine problem_reset_objectives(this)
676 do i = 1, this%n_objectives
677 call this%objective_list(i)%objective%reset_value()
679 end subroutine problem_reset_objectives
685 subroutine problem_reset_constraints(this)
689 do i = 1, this%n_constraints
690 call this%constraint_list(i)%constraint%reset_value()
692 end subroutine problem_reset_constraints
698 subroutine problem_reset_objective_sensitivities(this)
702 do i = 1, this%n_objectives
703 call this%objective_list(i)%objective%reset_sensitivity()
705 end subroutine problem_reset_objective_sensitivities
711 subroutine problem_reset_constraint_sensitivities(this)
715 do i = 1, this%n_constraints
716 call this%constraint_list(i)%constraint%reset_sensitivity()
718 end subroutine problem_reset_constraint_sensitivities
729 subroutine problem_accumulate_objectives(this, design, time)
731 class(
design_t),
intent(in) :: design
732 type(time_state_t),
intent(in) :: time
735 do i = 1, this%n_objectives
736 call this%objective_list(i)%objective%accumulate_value(
design, time)
738 end subroutine problem_accumulate_objectives
746 subroutine problem_accumulate_constraints(this, design, time)
748 class(
design_t),
intent(in) :: design
749 type(time_state_t),
intent(in) :: time
752 do i = 1, this%n_constraints
753 call this%constraint_list(i)%constraint%accumulate_value(
design, time)
755 end subroutine problem_accumulate_constraints
763 subroutine problem_accumulate_objective_sensitivities(this, design, time)
765 class(
design_t),
intent(in) :: design
766 type(time_state_t),
intent(in) :: time
769 do i = 1, this%n_objectives
770 call this%objective_list(i)%objective%accumulate_sensitivity(
design, &
773 end subroutine problem_accumulate_objective_sensitivities
781 subroutine problem_accumulate_constraint_sensitivities(this, design, time)
783 class(
design_t),
intent(in) :: design
784 type(time_state_t),
intent(in) :: time
787 do i = 1, this%n_constraints
788 call this%constraint_list(i)%constraint%accumulate_sensitivity(
design, &
791 end subroutine problem_accumulate_constraint_sensitivities
802 subroutine problem_get_objective_value(this, objective_value)
804 real(kind=rp),
intent(out) :: objective_value
807 objective_value = 0.0_rp
808 do i = 1, this%n_objectives
809 objective_value = objective_value + &
810 this%objective_list(i)%objective%get_weight() * &
811 this%objective_list(i)%objective%get_value()
814 end subroutine problem_get_objective_value
822 subroutine problem_get_all_objective_values(this, all_objective_values)
824 type(vector_t),
intent(inout) :: all_objective_values
827 do i = 1, this%n_objectives
828 all_objective_values%x(i) = this%objective_list(i)%objective%value
831 call all_objective_values%copy_from(host_to_device, sync = .true.)
833 end subroutine problem_get_all_objective_values
841 subroutine problem_get_constraint_values(this, constraint_value)
843 type(vector_t),
intent(inout) :: constraint_value
846 do i = 1, this%n_constraints
847 constraint_value%x(i) = this%constraint_list(i)%constraint%value
850 call constraint_value%copy_from(host_to_device, sync = .true.)
852 end subroutine problem_get_constraint_values
860 subroutine problem_get_objective_sensitivities(this, sensitivity)
862 type(vector_t),
intent(inout) :: sensitivity
865 call vector_cfill(sensitivity, 0.0_rp)
866 do i = 1, this%n_objectives
867 call vector_add2(sensitivity, &
868 this%objective_list(i)%objective%sensitivity)
871 end subroutine problem_get_objective_sensitivities
879 subroutine problem_get_constraint_sensitivities(this, sensitivity)
881 type(matrix_t),
target,
intent(inout) :: sensitivity
882 real(kind=rp),
pointer :: row(:)
886 do i = 1, this%n_constraints
887 call this%constraint_list(i)%constraint%sensitivity%copy_from( &
888 device_to_host, sync = i .eq. this%n_constraints)
891 do i = 1, this%n_constraints
892 row(1:this%n_design) => sensitivity%x(i, :)
894 call copy(row, this%constraint_list(i)%constraint%sensitivity%x, &
898 call sensitivity%copy_from(host_to_device, sync = .true.)
900 end subroutine problem_get_constraint_sensitivities
906 pure function problem_get_num_objectives(this)
result(n)
910 n = this%n_objectives
911 end function problem_get_num_objectives
914 pure function problem_get_num_constraints(this)
result(n)
918 n = this%n_constraints
919 end function problem_get_num_constraints
925 function problem_get_log_header(this, include_constraints)
result(buff)
927 logical,
intent(in),
optional :: include_constraints
928 character(len=4096) :: buff
929 character(len=128) :: mini_buff
930 character(len=128),
allocatable :: headers(:)
932 logical :: do_constraints
934 buff =
"Total objective function"
935 if (
present(include_constraints))
then
936 do_constraints = include_constraints
938 do_constraints = .true.
941 do i = 1, this%get_n_objectives()
942 n = this%objective_list(i)%objective%get_log_size()
945 call this%objective_list(i)%objective%get_log_headers(headers)
948 write(mini_buff,
'(", ", A)') trim(headers(j))
949 buff = trim(buff) // trim(mini_buff)
955 if (do_constraints)
then
956 do i = 1, this%get_n_constraints()
957 n = this%constraint_list(i)%constraint%get_log_size()
960 call this%constraint_list(i)%constraint%get_log_headers(headers)
963 write(mini_buff,
'(", ", A)') trim(headers(j))
964 buff = trim(buff) // trim(mini_buff)
971 end function problem_get_log_header
977 function problem_get_log_size(this, include_constraints)
result(n)
979 logical,
intent(in),
optional :: include_constraints
981 logical :: do_constraints
984 if (
present(include_constraints))
then
985 do_constraints = include_constraints
987 do_constraints = .true.
989 do i = 1, this%get_n_objectives()
990 n = n + this%objective_list(i)%objective%get_log_size()
993 if (do_constraints)
then
994 do i = 1, this%get_n_constraints()
995 n = n + this%constraint_list(i)%constraint%get_log_size()
999 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.