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, log_size
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
111 procedure, pass(this) :: check_objective_windows => &
112 problem_check_objective_windows
117 procedure, pass(this),
public :: read_objectives => problem_read_objectives
119 procedure, pass(this),
public :: read_constraints => &
120 problem_read_constraints
126 procedure, pass(this),
public :: write => problem_write
128 procedure, pass(this),
public :: get_log_values => problem_get_log_values
131 procedure, pass(this),
public :: add_objective => problem_add_objective
133 procedure, pass(this),
public :: add_constraint => problem_add_constraint
139 procedure, pass(this) :: update_objectives => &
140 problem_update_objectives
142 procedure, pass(this) :: update_constraints => &
143 problem_update_constraints
145 procedure, pass(this) :: update_objective_sensitivities => &
146 problem_update_objective_sensitivities
148 procedure, pass(this) :: update_constraint_sensitivities => &
149 problem_update_constraint_sensitivities
152 procedure, pass(this) :: reset_objectives => &
153 problem_reset_objectives
155 procedure, pass(this) :: reset_constraints => &
156 problem_reset_constraints
158 procedure, pass(this) :: reset_objective_sensitivities => &
159 problem_reset_objective_sensitivities
161 procedure, pass(this) :: reset_constraint_sensitivities => &
162 problem_reset_constraint_sensitivities
165 procedure, pass(this) :: accumulate_objectives => &
166 problem_accumulate_objectives
168 procedure, pass(this) :: accumulate_constraints => &
169 problem_accumulate_constraints
171 procedure, pass(this) :: accumulate_objective_sensitivities => &
172 problem_accumulate_objective_sensitivities
174 procedure, pass(this) :: accumulate_constraint_sensitivities => &
175 problem_accumulate_constraint_sensitivities
181 procedure, pass(this),
public :: get_objective_value => &
182 problem_get_objective_value
184 procedure, pass(this),
public :: get_all_objective_values => &
185 problem_get_all_objective_values
187 procedure, pass(this),
public :: get_constraint_values => &
188 problem_get_constraint_values
190 procedure, pass(this),
public :: get_objective_sensitivities => &
191 problem_get_objective_sensitivities
193 procedure, pass(this),
public :: get_constraint_sensitivities => &
194 problem_get_constraint_sensitivities
197 procedure, pass(this) :: get_n_objectives => problem_get_num_objectives
199 procedure, pass(this) :: get_n_constraints => problem_get_num_constraints
202 procedure, pass(this) :: get_log_header => problem_get_log_header
204 procedure, pass(this) :: get_log_size => problem_get_log_size
216 type(json_file),
intent(inout) :: parameters
217 class(
design_t),
intent(in) :: design
218 type(
simulation_t),
optional,
intent(inout) :: simulation
222 this%n_design =
design%size()
225 call this%read_objectives(parameters,
design, simulation)
226 call this%read_constraints(parameters,
design, simulation)
231 subroutine problem_free(this)
236 this%n_objectives = 0
237 this%n_constraints = 0
240 if (
allocated(this%objective_list))
then
241 do i = 1,
size(this%objective_list)
242 call this%objective_list(i)%free()
244 deallocate(this%objective_list)
248 if (
allocated(this%constraint_list))
then
249 do i = 1,
size(this%constraint_list)
250 call this%constraint_list(i)%free()
252 deallocate(this%constraint_list)
254 end subroutine problem_free
257 subroutine problem_write(this, idx)
259 integer,
intent(in) :: idx
261 end subroutine problem_write
267 subroutine problem_get_log_values(this, values, include_constraints)
269 real(kind=rp),
intent(out) :: values(:)
270 logical,
intent(in),
optional :: include_constraints
271 integer :: i, n, offset
272 real(kind=rp) :: objective_value
273 real(kind=rp),
allocatable :: tmp(:)
274 logical :: do_constraints
276 if (
present(include_constraints))
then
277 do_constraints = include_constraints
279 do_constraints = .true.
282 call this%get_objective_value(objective_value)
284 values(1) = objective_value
287 do i = 1, this%n_objectives
288 n = this%objective_list(i)%objective%get_log_size()
291 call this%objective_list(i)%objective%get_log_values(tmp)
292 values(offset:offset + n - 1) = tmp
298 if (do_constraints)
then
299 do i = 1, this%n_constraints
300 n = this%constraint_list(i)%constraint%get_log_size()
303 call this%constraint_list(i)%constraint%get_log_values(tmp)
304 values(offset:offset + n - 1) = tmp
310 end subroutine problem_get_log_values
316 subroutine problem_read_objectives(this, parameters, design, simulation)
318 type(json_file),
intent(inout) :: parameters
319 class(
design_t),
intent(in) :: design
320 type(
simulation_t),
optional,
intent(inout) :: simulation
324 character(len=:),
allocatable :: path, type
325 type(json_file) :: objective_json
326 integer :: n_objectives, i
329 call neko_log%section(
"Reading objectives")
332 path =
"optimization.objectives"
333 if (parameters%valid_path(path))
then
334 call parameters%info(path, n_children = n_objectives)
337 do i = 1, n_objectives
338 call json_extract_item(parameters, path, i, objective_json)
339 call json_get(objective_json,
"type", type)
340 call neko_log%message(type)
347 if (
present(simulation))
then
352 call json_get_or_default(parameters, &
353 "adjoint_fluid.dealias_sensitivity", dealias, .true.)
354 call alo%init_from_attributes(
design, simulation, weight = 1.0_rp, &
355 name =
"Augmented Lagrangian", mask_name =
"", &
361 call neko_log%end_section()
363 end subroutine problem_read_objectives
366 subroutine problem_read_constraints(this, parameters, design, simulation)
368 type(json_file),
intent(inout) :: parameters
369 class(
design_t),
intent(in) :: design
371 type(
simulation_t),
optional,
intent(inout) :: simulation
374 character(len=:),
allocatable :: path, type
375 type(json_file) :: constraint_json
376 integer :: n_constraints, i
378 call neko_log%section(
"Reading constraints")
381 path =
"optimization.constraints"
383 if (parameters%valid_path(path))
then
384 call parameters%info(path, n_children = n_constraints)
387 do i = 1, n_constraints
388 call json_extract_item(parameters, path, i, constraint_json)
389 call json_get(constraint_json,
"type", type)
390 call neko_log%message(type)
398 call neko_log%end_section()
400 end subroutine problem_read_constraints
403 subroutine problem_add_objective(this, objective)
405 class(
objective_t),
allocatable,
intent(inout) :: objective
410 if (
allocated(this%objective_list))
then
411 n =
size(this%objective_list)
412 call move_alloc(this%objective_list, temp_list)
413 allocate(this%objective_list(n + 1))
414 if (
allocated(temp_list))
then
416 call move_alloc(temp_list(i)%objective, &
417 this%objective_list(i)%objective)
421 allocate(this%objective_list(1))
424 call move_alloc(
objective, this%objective_list(n + 1)%objective)
425 this%n_objectives = n + 1
426 end subroutine problem_add_objective
429 subroutine problem_add_constraint(this, constraint)
431 class(
constraint_t),
allocatable,
intent(inout) :: constraint
436 if (
allocated(this%constraint_list))
then
437 n =
size(this%constraint_list)
438 call move_alloc(this%constraint_list, temp_list)
439 allocate(this%constraint_list(n + 1))
440 if (
allocated(temp_list))
then
442 call move_alloc(temp_list(i)%constraint, &
443 this%constraint_list(i)%constraint)
447 allocate(this%constraint_list(1))
450 call move_alloc(
constraint, this%constraint_list(n + 1)%constraint)
451 this%n_constraints = n + 1
452 end subroutine problem_add_constraint
458 subroutine problem_compute(this, design, simulation)
460 class(
design_t),
intent(inout) :: design
461 class(
simulation_t),
optional,
intent(inout) :: simulation
463 if (
present(simulation))
then
464 call simulation%reset()
465 if (simulation%unsteady)
then
467 call this%run_forward_unsteady(simulation,
design)
469 call simulation%run_forward()
471 call this%update_objectives(
design)
474 call this%update_objectives(
design)
477 call this%update_constraints(
design)
479 end subroutine problem_compute
482 subroutine problem_compute_sensitivity(this, design, simulation)
484 class(
design_t),
intent(inout) :: design
485 class(
simulation_t),
optional,
intent(inout) :: simulation
487 type(vector_t) :: objective_sensitivity
489 if (
present(simulation))
then
490 if (simulation%unsteady)
then
492 call this%run_backward_unsteady(simulation,
design)
494 call simulation%run_backward()
496 call this%update_objective_sensitivities(
design)
499 call this%update_objective_sensitivities(
design)
502 call this%update_constraint_sensitivities(
design)
504 call objective_sensitivity%init(this%n_design)
505 call this%get_objective_sensitivities(objective_sensitivity)
507 call design%map_backward(objective_sensitivity)
509 call objective_sensitivity%free()
510 end subroutine problem_compute_sensitivity
513 subroutine problem_run_forward_unsteady(this, simulation, design)
516 class(
design_t),
intent(inout) :: design
517 type(time_step_controller_t) :: dt_controller
518 real(kind=dp) :: loop_start
520 call dt_controller%init(simulation%neko_case%params)
522 call simulation%reset()
523 call simulation_init(simulation%neko_case, dt_controller)
526 call this%reset_objectives()
528 if (.not.
allocated(simulation%state_recover))
then
529 call neko_error(
"State recovery not initialized.")
532 call profiler_start_region(
"Forward simulation")
533 loop_start = mpi_wtime()
534 simulation%n_timesteps = 0
535 do while (simulation%neko_case%time%t .lt. &
536 simulation%neko_case%time%end_time)
537 simulation%n_timesteps = simulation%n_timesteps + 1
539 call simulation_step(simulation%neko_case, dt_controller, loop_start)
541 call this%accumulate_objectives(
design, simulation%neko_case%time)
543 call simulation%state_recover%save()
545 call profiler_end_region(
"Forward simulation")
547 call this%check_objective_windows()
549 call simulation_finalize(simulation%neko_case)
551 end subroutine problem_run_forward_unsteady
558 subroutine problem_check_objective_windows(this)
560 character(len=LOG_SIZE) :: log_buf
563 do i = 1, this%n_objectives
564 if (this%objective_list(i)%objective%value_weight .gt. 0.0_rp) cycle
566 write (log_buf,
'(A,A,A)')
"Objective '", &
567 trim(this%objective_list(i)%objective%name), &
568 "' was never sampled; its time window misses the run."
569 call neko_log%warning(trim(log_buf))
571 end subroutine problem_check_objective_windows
574 subroutine problem_run_backward_unsteady(this, simulation, design)
577 class(
design_t),
intent(inout) :: design
578 type(time_step_controller_t) :: dt_controller
579 real(kind=dp) :: loop_start
581 real(kind=rp) :: total_time
583 type(time_state_t) :: accumulation_time
585 call dt_controller%init(simulation%neko_case%params)
590 call this%reset_objective_sensitivities()
592 cfl = simulation%adjoint_case%fluid_adj%compute_cfl( &
593 simulation%adjoint_case%time%dt)
594 loop_start = mpi_wtime()
596 if (.not.
allocated(simulation%state_recover))
then
597 call neko_error(
"State recovery not initialized.")
601 total_time = simulation%n_timesteps * simulation%adjoint_case%time%dt
603 call profiler_start_region(
"Adjoint simulation")
605 do i = simulation%n_timesteps, 1, -1
607 call simulation%state_recover%restore(i)
609 accumulation_time = simulation%adjoint_case%time
610 accumulation_time%t = total_time - simulation%adjoint_case%time%t
611 call this%accumulate_objective_sensitivities(
design, accumulation_time)
614 cfl, loop_start, total_time)
617 call profiler_end_region(
"Adjoint simulation")
621 end subroutine problem_run_backward_unsteady
632 subroutine problem_update_objectives(this, design)
634 class(
design_t),
intent(in) :: design
637 do i = 1, this%n_objectives
638 call this%objective_list(i)%objective%update_value(
design)
640 end subroutine problem_update_objectives
648 subroutine problem_update_constraints(this, design)
650 class(
design_t),
intent(in) :: design
653 do i = 1, this%n_constraints
654 call this%constraint_list(i)%constraint%update_value(
design)
656 end subroutine problem_update_constraints
664 subroutine problem_update_objective_sensitivities(this, design)
666 class(
design_t),
intent(in) :: design
669 do i = 1, this%n_objectives
670 call this%objective_list(i)%objective%update_sensitivity(
design)
672 end subroutine problem_update_objective_sensitivities
680 subroutine problem_update_constraint_sensitivities(this, design)
682 class(
design_t),
intent(in) :: design
685 do i = 1, this%n_constraints
686 call this%constraint_list(i)%constraint%update_sensitivity(
design)
688 end subroutine problem_update_constraint_sensitivities
697 subroutine problem_reset_objectives(this)
701 do i = 1, this%n_objectives
702 call this%objective_list(i)%objective%reset_value()
704 end subroutine problem_reset_objectives
710 subroutine problem_reset_constraints(this)
714 do i = 1, this%n_constraints
715 call this%constraint_list(i)%constraint%reset_value()
717 end subroutine problem_reset_constraints
723 subroutine problem_reset_objective_sensitivities(this)
727 do i = 1, this%n_objectives
728 call this%objective_list(i)%objective%reset_sensitivity()
730 end subroutine problem_reset_objective_sensitivities
736 subroutine problem_reset_constraint_sensitivities(this)
740 do i = 1, this%n_constraints
741 call this%constraint_list(i)%constraint%reset_sensitivity()
743 end subroutine problem_reset_constraint_sensitivities
754 subroutine problem_accumulate_objectives(this, design, time)
756 class(
design_t),
intent(in) :: design
757 type(time_state_t),
intent(in) :: time
760 do i = 1, this%n_objectives
761 call this%objective_list(i)%objective%accumulate_value(
design, time)
763 end subroutine problem_accumulate_objectives
771 subroutine problem_accumulate_constraints(this, design, time)
773 class(
design_t),
intent(in) :: design
774 type(time_state_t),
intent(in) :: time
777 do i = 1, this%n_constraints
778 call this%constraint_list(i)%constraint%accumulate_value(
design, time)
780 end subroutine problem_accumulate_constraints
788 subroutine problem_accumulate_objective_sensitivities(this, design, time)
790 class(
design_t),
intent(in) :: design
791 type(time_state_t),
intent(in) :: time
794 do i = 1, this%n_objectives
795 call this%objective_list(i)%objective%accumulate_sensitivity(
design, &
798 end subroutine problem_accumulate_objective_sensitivities
806 subroutine problem_accumulate_constraint_sensitivities(this, design, time)
808 class(
design_t),
intent(in) :: design
809 type(time_state_t),
intent(in) :: time
812 do i = 1, this%n_constraints
813 call this%constraint_list(i)%constraint%accumulate_sensitivity(
design, &
816 end subroutine problem_accumulate_constraint_sensitivities
827 subroutine problem_get_objective_value(this, objective_value)
829 real(kind=rp),
intent(out) :: objective_value
832 objective_value = 0.0_rp
833 do i = 1, this%n_objectives
834 objective_value = objective_value + &
835 this%objective_list(i)%objective%get_weight() * &
836 this%objective_list(i)%objective%get_value()
839 end subroutine problem_get_objective_value
847 subroutine problem_get_all_objective_values(this, all_objective_values)
849 type(vector_t),
intent(inout) :: all_objective_values
852 do i = 1, this%n_objectives
853 all_objective_values%x(i) = this%objective_list(i)%objective%value
856 call all_objective_values%copy_from(host_to_device, sync = .true.)
858 end subroutine problem_get_all_objective_values
866 subroutine problem_get_constraint_values(this, constraint_value)
868 type(vector_t),
intent(inout) :: constraint_value
871 do i = 1, this%n_constraints
872 constraint_value%x(i) = this%constraint_list(i)%constraint%value
875 call constraint_value%copy_from(host_to_device, sync = .true.)
877 end subroutine problem_get_constraint_values
885 subroutine problem_get_objective_sensitivities(this, sensitivity)
887 type(vector_t),
intent(inout) :: sensitivity
890 call vector_cfill(sensitivity, 0.0_rp)
891 do i = 1, this%n_objectives
892 call vector_add2(sensitivity, &
893 this%objective_list(i)%objective%sensitivity)
896 end subroutine problem_get_objective_sensitivities
904 subroutine problem_get_constraint_sensitivities(this, sensitivity)
906 type(matrix_t),
target,
intent(inout) :: sensitivity
907 real(kind=rp),
pointer :: row(:)
911 do i = 1, this%n_constraints
912 call this%constraint_list(i)%constraint%sensitivity%copy_from( &
913 device_to_host, sync = i .eq. this%n_constraints)
916 do i = 1, this%n_constraints
917 row(1:this%n_design) => sensitivity%x(i, :)
919 call copy(row, this%constraint_list(i)%constraint%sensitivity%x, &
923 call sensitivity%copy_from(host_to_device, sync = .true.)
925 end subroutine problem_get_constraint_sensitivities
931 pure function problem_get_num_objectives(this)
result(n)
935 n = this%n_objectives
936 end function problem_get_num_objectives
939 pure function problem_get_num_constraints(this)
result(n)
943 n = this%n_constraints
944 end function problem_get_num_constraints
950 function problem_get_log_header(this, include_constraints)
result(buff)
952 logical,
intent(in),
optional :: include_constraints
953 character(len=4096) :: buff
954 character(len=128) :: mini_buff
955 character(len=128),
allocatable :: headers(:)
957 logical :: do_constraints
959 buff =
"Total objective function"
960 if (
present(include_constraints))
then
961 do_constraints = include_constraints
963 do_constraints = .true.
966 do i = 1, this%get_n_objectives()
967 n = this%objective_list(i)%objective%get_log_size()
970 call this%objective_list(i)%objective%get_log_headers(headers)
973 write(mini_buff,
'(", ", A)') trim(headers(j))
974 buff = trim(buff) // trim(mini_buff)
980 if (do_constraints)
then
981 do i = 1, this%get_n_constraints()
982 n = this%constraint_list(i)%constraint%get_log_size()
985 call this%constraint_list(i)%constraint%get_log_headers(headers)
988 write(mini_buff,
'(", ", A)') trim(headers(j))
989 buff = trim(buff) // trim(mini_buff)
996 end function problem_get_log_header
1002 function problem_get_log_size(this, include_constraints)
result(n)
1004 logical,
intent(in),
optional :: include_constraints
1006 logical :: do_constraints
1009 if (
present(include_constraints))
then
1010 do_constraints = include_constraints
1012 do_constraints = .true.
1014 do i = 1, this%get_n_objectives()
1015 n = n + this%objective_list(i)%objective%get_log_size()
1018 if (do_constraints)
then
1019 do i = 1, this%get_n_constraints()
1020 n = n + this%constraint_list(i)%constraint%get_log_size()
1024 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.