← Back to branch summary

~fluidity-core/fluidity/adjoint

« back to all changes in this revision

Viewing changes to adjoint/Adjoint_Controls.F90

Merge the fluidity_revolve branch.

Show diffs side-by-side

added added

removed removed

Lines of Context:
adjoint/Adjoint_Controls.F90
43
43
 
44
44
    contains 
45
45
 
46
 
    ! Retrieves the control details from control number control_no in the optin tree. The outputs are:
 
46
 
 
47
    ! Returns the material phase and field name of the control_no'th control
 
48
    subroutine get_control_field_name(control_no, material_phase_name, field_name)
 
49
      integer, intent(in) :: control_no
 
50
      character(len=OPTION_PATH_LEN), intent(out) :: field_name, material_phase_name
 
51
 
 
52
      character(len=OPTION_PATH_LEN) :: control_type, control_name, name 
 
53
      integer :: s_idx
 
54
 
 
55
      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/name", control_name)
 
56
      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/type/name", control_type)
 
57
      call get_option("/adjoint/controls/control[" // int2str(control_no) //"]/type::" // trim(control_type) // "/field_name", name)
 
58
 
 
59
      s_idx = scan(trim(name), "::")
 
60
      if (s_idx == 0) then 
 
61
        FLExit("The control " // trim(control_name) // " uses an invalid field_name. It should be of the form Materialphase::Field")
 
62
      end if
 
63
      material_phase_name = name(1:s_idx - 1)
 
64
      field_name = name(s_idx + 2:len_trim(name))
 
65
    end subroutine get_control_field_name
 
66
 
 
67
    ! Retrieves the control details from control number control_no in the option tree. The outputs are:
47
68
    subroutine get_control_details(states, timestep, control_no, control_name, state_id, field_name, field_type, &
48
69
                                 & active, have_lb, lb_state_id, lb_material_phase_name, lb_field_name, have_ub, &
49
70
                                 & ub_state_id, ub_material_phase_name, ub_field_name)
58
79
                                                                                          ! field_name: the name of the associated field
59
80
                                                                                          ! field_type: the type of the associated field
60
81
      integer, intent(out) :: state_id ! the state number in which the associated control field exists
61
 
      logical, intent(out) :: active ! false if the control is active in this timestep (e.g. InitialConditions for timesteps >0)
 
82
      logical, intent(out) :: active ! false if the control is not active in this timestep (e.g. InitialConditions for timesteps >0)
62
83
      logical, intent(out), optional :: have_ub, have_lb
63
84
      integer, intent(out), optional :: lb_state_id, ub_state_id
64
85
      character(len=OPTION_PATH_LEN), intent(out), optional :: lb_material_phase_name, ub_material_phase_name
89
110
          exit
90
111
        end if
91
112
      end do
92
 
      if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
 
113
      if (state_id==size(states)+1) then
93
114
        FLExit("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_name) // ".")
94
115
      end if
95
116
      
163
184
          lb_field_name = name(s_idx + 2:len_trim(name))
164
185
          ! Find state associated with the lower bound material phase
165
186
          do lb_state_id = 1, size(states) 
166
 
            if (trim(states(lb_state_id)%name) == trim(lb_material_phase_name)) then
 
187
            if (trim(states(lb_state_id)%name) .eq. trim(lb_material_phase_name)) then
167
188
              exit
168
189
            end if
169
190
          end do
170
 
          if (.not. trim(states(lb_state_id)%name) == trim(lb_material_phase_name)) then
 
191
          if (.not. trim(states(lb_state_id)%name) .eq. trim(lb_material_phase_name)) then
171
192
            FLExit("Could not find state " // trim(lb_material_phase_name) // " as specified in the lower bound of control " // trim(control_name) // ".")
172
193
          end if
173
194
        end if
368
389
    end subroutine adjoint_write_controls
369
390
 
370
391
    ! The counterpart of adjoint_write_controls: Loop through all controls in the option file and fill the fields from the control files.
371
 
    subroutine adjoint_load_controls(timestep, dt, states)
 
392
    ! If a field name filter is provided, then only the field with the specified name will be loaded
 
393
    subroutine adjoint_load_controls(timestep, dt, states, field_name_filter)
372
394
      integer, intent(in) :: timestep
373
395
      real, intent(in) :: dt
374
396
      type(state_type), dimension(:), intent(in) :: states
 
397
      character(len=*), intent(in), optional :: field_name_filter 
375
398
#ifdef HAVE_ADJOINT
376
399
      integer :: nb_controls, i, state_id, s_idx
377
 
      logical :: active, file_exists
378
 
      character(len=OPTION_PATH_LEN) :: field_name, field_type, control_name, simulation_name, control_type
 
400
      logical :: active, file_exists, filter_active
 
401
      character(len=OPTION_PATH_LEN) :: material_phase_name, field_name, field_type, control_name, simulation_name, control_type
379
402
 
380
403
      ! Do not read anything if we are supposed to write the controls
381
404
      if (.not. have_option("/adjoint/controls/load_controls")) then
387
410
 
388
411
      ! Now loop over the controls, read their controls files and fill the associated fields
389
412
      do i = 0, nb_controls-1
 
413
        call get_control_field_name(i, material_phase_name, field_name)
 
414
        if (present(field_name_filter)) then
 
415
          if (trim(field_name) /= trim(field_name_filter)) then
 
416
            cycle 
 
417
          end if
 
418
        end if
 
419
 
390
420
        call get_control_details(states, timestep, i, control_name, state_id, field_name, field_type, active)
 
421
 
391
422
        call get_option("/adjoint/controls/control[" // int2str(i) // "]/type/name", control_type)
392
423
 
393
424
        ! We load all OTHER controls, then apply our strong boundary conditions, then load BC controls
394
425
        if (trim(control_type) == "boundary_condition") cycle
395
426
 
396
427
        if (active) then
397
 
          ! Check that the control parameter file extists
 
428
          ! Check that the control parameter file exists
398
429
          inquire(file="control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl", exist=file_exists)
399
430
          if (.not. file_exists) then
400
431
            ewrite (0, *) "File 'control_" // trim(simulation_name) // "_" // trim(control_name) // &
418
449
      call set_dirichlet_consistent(states)
419
450
 
420
451
      do i = 0, nb_controls-1
 
452
        call get_control_field_name(i, material_phase_name, field_name)
 
453
        if (present(field_name_filter)) then
 
454
          if (trim(field_name) /= trim(field_name_filter)) then
 
455
            cycle 
 
456
          end if
 
457
        end if
 
458
 
421
459
        call get_control_details(states, timestep, i, control_name, state_id, field_name, field_type, active)
 
460
 
422
461
        call get_option("/adjoint/controls/control[" // int2str(i) // "]/type/name", control_type)
423
462
 
424
463
        ! We load all OTHER controls, then apply our strong boundary conditions, then load BC controls
425
464
        if (trim(control_type) /= "boundary_condition") cycle
426
465
 
427
466
        if (active) then
428
 
          ! Check that the control parameter file extists
 
467
          ! Check that the control parameter file exists
429
468
          inquire(file="control_" // trim(simulation_name) // "_" // trim(control_name) // "_" // int2str(timestep) // ".pkl", exist=file_exists)
430
469
          if (.not. file_exists) then
431
470
            ewrite (0, *) "File 'control_" // trim(simulation_name) // "_" // trim(control_name) // &
466
505
      ! Remove the suffix _adjoint from the simulation_name
467
506
 
468
507
      nb_controls = option_count("/adjoint/controls/control")
 
508
 
469
509
      do i = 0, nb_controls-1
 
510
 
470
511
        call get_option("/adjoint/controls/control[" // int2str(i) //"]/name", control_deriv_name)
471
512
        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type/name", control_type)
472
513
        call get_option("/adjoint/controls/control[" // int2str(i) //"]/type::" // trim(control_type) // "/field_name", name)
 
514
 
473
515
        s_idx = scan(trim(name), "::")
474
516
        if (s_idx == 0) then 
475
517
          FLExit("The control " // trim(control_deriv_name) // " uses an invalid field_name. It should be of the form Materialphase::Fieldname")
476
518
        end if
477
519
        material_phase_name = name(1:s_idx - 1)
478
520
        field_name = name(s_idx + 2:len_trim(name))
 
521
 
479
522
        ! Find state associated with the material phase
480
523
        do state_id = 1, size(states) 
481
524
          if (trim(states(state_id)%name) == trim(material_phase_name)) then
482
525
            exit
483
526
          end if
484
527
        end do
 
528
 
485
529
        ! Make sure we found the state
486
 
        if (.not. trim(states(state_id)%name) == trim(material_phase_name)) then
 
530
        if (state_id == size(states)+1) then
487
531
          FLExit("Could not find state " // trim(material_phase_name) // " as specified in control " // trim(control_deriv_name) // ".")
488
532
        end if
489
533
 
515
559
        elseif (trim(control_type) == "boundary_condition") then
516
560
          FLAbort("Boundary condition control not implemented yet.")
517
561
        end if
 
562
 
518
563
        ! Save the control parameter to disk
519
564
        call python_reset()
520
565
        call python_add_state(states(state_id))