mma_simcomp.f90

Go to the documentation of this file.

! Copyright (c) 2024, The Neko Authors
 ! All rights reserved.
 !
 ! Redistribution and use in source and binary forms, with or without
 ! modification, are permitted provided that the following conditions
 ! are met:
 !
 !   * Redistributions of source code must retain the above copyright
 !     notice, this list of conditions and the following disclaimer.
 !
 !   * Redistributions in binary form must reproduce the above
 !     copyright notice, this list of conditions and the following
 !     disclaimer in the documentation and/or other materials provided
 !     with the distribution.
 !
 !   * Neither the name of the authors nor the names of its
 !     contributors may be used to endorse or promote products derived
 !     from this software without specific prior written permission.
 !
 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 ! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
 ! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
 ! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
 ! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 ! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 ! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
 ! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 ! POSSIBILITY OF SUCH DAMAGE.
 
 ! Implements the `mma_comp_t` type.
module mma_simcomp
  use num_types, only: rp
  use case, only: case_t
  use json_module, only: json_file
  use json_utils, only: json_get_or_default
  use simulation_component, only: simulation_component_t
  use field, only: field_t
  use logger, only: neko_log
  use vector, only: vector_t
  use matrix, only: matrix_t
  use mma, only: mma_t
 
  use device, only: device_memcpy, host_to_device, device_to_host
  use mpi_f08, only: mpi_allreduce, mpi_integer, mpi_sum, mpi_comm_world, &
       mpi_min, mpi_max, mpi_in_place
  use comm, only: pe_rank, neko_comm, mpi_real_precision
  implicit none
  private
 
  ! An empty user defined simulation component.
  ! This is a simple example of a user-defined simulation component.
  type, public, extends(simulation_component_t) :: mma_comp_t
 
     type(field_t) :: designx
 
     integer :: m
 
     real(kind=rp) :: scale
     logical :: auto_scale
 
     type(mma_t) :: mma
 
   contains
     ! Constructor from json, wrapping the actual constructor.
     procedure, pass(this) :: init => simcomp_test_init_from_json
     ! Destructor.
     procedure, pass(this) :: free => simcomp_test_free
     ! Compute the simcomp_test field.
     procedure, pass(this) :: compute_ => simcomp_test_compute
  end type mma_comp_t
 
contains
  ! Constructor from json.
  subroutine simcomp_test_init_from_json(this, json, case)
    class(mma_comp_t), intent(inout) :: this
    type(json_file), intent(inout) :: json
    class(case_t), intent(inout), target :: case
    integer :: nloc
 
    call this%init_base(json, case)
    call this%designx%init(case%msh, case%fluid%Xh, "designx")
 
    call json_get_or_default(json, "m", this%m, 2)
    nloc = this%designx%dof%size()
 
    ! initial design
    this%designx%x = 1.0
 
    call this%mma%init( reshape(this%designx%x, [nloc]), &
         nloc, this%m, case%params, this%scale, this%auto_scale)
  end subroutine simcomp_test_init_from_json
 
 
  ! Destructor.
  subroutine simcomp_test_free(this)
    class(mma_comp_t), intent(inout) :: this
    call this%mma%free()
    call this%free_base()
  end subroutine simcomp_test_free
 
  ! Computations.
  subroutine simcomp_test_compute(this, t, tstep)
    implicit none
    class(mma_comp_t), intent(inout) :: this
    real(kind=rp), intent(in) :: t
    integer, intent(in) :: tstep
 
    integer :: iter, i, ierr, rank, size, nglobal
    integer, allocatable :: recv_counts(:), displs(:)
 
    real(kind=rp) :: start_time, end_time
    real(kind=rp), dimension(this%mma%get_n()) :: x
    real(kind=rp) :: f0val
    type(vector_t) :: df0dx, fval
    type(matrix_t) :: dfdx
 
    real(kind=rp), dimension(this%mma%get_n(), 4) :: stuff
 
 
    real(kind=rp), allocatable :: all_stuff(:, :)
    integer, allocatable :: nloc_all(:)
 
    call df0dx%init(this%mma%get_n())
    call fval%init(this%mma%get_m())
    call dfdx%init(this%mma%get_m(), this%mma%get_n())
 
 
    call cpu_time(start_time)
    x = reshape(this%designx%x, [this%mma%get_n()])
    call func1(this, this%mma%get_n(), this%mma%get_m(), &
         f0val, df0dx%x, fval%x, dfdx%x)
    ! update the device pointer
    call device_memcpy(df0dx%x, df0dx%x_d, this%mma%get_n(), &
         host_to_device, sync = .false.)
    call device_memcpy(fval%x, fval%x_d, this%mma%get_m(), &
         host_to_device, sync = .false.)
    call device_memcpy(dfdx%x, dfdx%x_d, this%mma%get_n()*this%mma%get_m(), &
         host_to_device, sync = .false.)
 
    if (pe_rank .eq. 0) then
       print *, 'iter = ', 0, &
            '-------, f0val = ', f0val, ',   fval = ', fval%x
    end if
 
    stuff(:, 1) = reshape(this%designx%dof%x, [this%mma%get_n()])
    stuff(:, 2) = reshape(this%designx%dof%y, [this%mma%get_n()])
    stuff(:, 3) = reshape(this%designx%dof%z, [this%mma%get_n()])
    stuff(:, 4) = reshape(this%designx%x, [this%mma%get_n()])
 
    call mpi_comm_size(neko_comm, size, ierr)
    call mpi_comm_rank(neko_comm, rank, ierr)
    allocate(recv_counts(size))
    allocate(displs(size))
    call mpi_allreduce(this%mma%get_n(), nglobal, 1, &
         mpi_integer, mpi_sum, neko_comm, ierr)
 
    allocate(nloc_all(size))
    !!!!Use MPI_Allgather to gather the `nloc` from each process into `nloc_all`
    call mpi_allgather(this%mma%get_n(), 1, mpi_integer, nloc_all, &
         1, mpi_integer, mpi_comm_world, ierr)
    recv_counts = this%mma%get_n()
    displs(1) = 0
    do i = 2, size
       displs(i) = displs(i-1) + nloc_all(i-1)
    end do
    allocate(all_stuff(nglobal, 4))
    ! Gather data from all processes to the root process
    call mpi_gatherv(stuff(:, 1), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 1), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 2), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 2), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 3), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 3), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 4), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 4), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    ! Only root process writes the file
    if (pe_rank .eq. 0) then
       call write_stuff_vtk(all_stuff, nglobal, "stuff_init.vtk")
    end if
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! The optimization loop
    do iter = 1, 100 !10
       call this%mma%update(iter, x, df0dx, fval, dfdx)
       this%designx%x = reshape(x, shape(this%designx%x))
 
       call func1(this, this%mma%get_n(), this%mma%get_m(), &
            f0val, df0dx%x, fval%x, dfdx%x)
       ! update the device pointer
       call device_memcpy(df0dx%x, df0dx%x_d, this%mma%get_n(), &
            host_to_device, sync = .false.)
       call device_memcpy(fval%x, fval%x_d, this%mma%get_m(), &
            host_to_device, sync = .false.)
       call device_memcpy(dfdx%x, dfdx%x_d, &
            this%mma%get_n()*this%mma%get_m(), host_to_device, sync = .false.)
 
       call this%mma%KKT(this%designx%x, df0dx, fval, dfdx)
 
       if (pe_rank .eq. 0) then
          print *, 'iter = ', iter, &
               '-------, f0val = ', f0val, ',   fval = ', fval%x(1), &
               ',  KKTmax = ', this%mma%get_residumax(), &
               ', KKTnorm2 = ', this%mma%get_residunorm()
       end if
 
       if (this%mma%get_residumax() .lt. 1.0e-3_rp) exit
    end do
 
    call cpu_time(end_time)
 
    print *, 'Elapsed Time: ', end_time - start_time, ' seconds'
    print *, "this%designx%x is updated"
 
 
    stuff(:, 1) = reshape(this%designx%dof%x, [this%mma%get_n()])
    stuff(:, 2) = reshape(this%designx%dof%y, [this%mma%get_n()])
    stuff(:, 3) = reshape(this%designx%dof%z, [this%mma%get_n()])
    stuff(:, 4) = reshape(this%designx%x, [this%mma%get_n()])
 
    call mpi_gatherv(stuff(:, 1), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 1), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 2), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 2), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 3), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 3), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    call mpi_gatherv(stuff(:, 4), this%mma%get_n(), mpi_real_precision, &
         all_stuff(:, 4), recv_counts, displs, mpi_real_precision, &
         0, neko_comm, ierr)
    ! Only root process writes the file
    if (pe_rank .eq. 0) then
       call write_stuff_vtk(all_stuff, nglobal, "stuff_final.vtk")
    end if
 
  end subroutine simcomp_test_compute
 
  subroutine write_stuff_vtk(stuff, n, filename)
    ! ----------------------------------------------------------- !
    !  This subroutine writes a nx4 array into a vtk file.        !
    !  The array is called stuff(n, 4) holding the coordinates of !
    !  all points and thier corresponding scalar value            !
    !  For a given point n, the array is defined as follows:      !
    !  x = stuff(n, 1), y = stuff(n, 2), z = stuff(n, 3)          !
    !  scalar field value = stuff(n, 4)                           !
    ! ----------------------------------------------------------- !
    integer, intent(in) :: n ! Number of design variables
    real(kind=rp), dimension(n, 4), intent(in) :: stuff
    character(len=*), intent(in) :: filename ! Filename to save the VTK data
 
    integer :: i ! Loop variable
    integer :: unit ! File unit number
 
    ! Open the VTK file for writing
    open(newunit = unit, file = filename, status = "replace", action = "write")
 
    ! Write VTK header
    write(unit, '(A)') '# vtk DataFile Version 3.0'
    write(unit, '(A)') 'Stuff data'
    write(unit, '(A)') 'ASCII'
    write(unit, '(A)') 'DATASET POLYDATA'
 
    ! Write the points (x, y, z coordinates)
    write(unit, '(A, I8)') 'POINTS ', n, ' double'
    do i = 1, n
       write(unit, '(3(ES15.8, 1X))') stuff(i, 1), stuff(i, 2), stuff(i, 3)
    end do
 
    ! Write the temperature data associated with the points
    write(unit, '(A, I8)') 'POINT_DATA ', n
    write(unit, '(A)') 'SCALARS Temperature double 1'
    write(unit, '(A)') 'LOOKUP_TABLE default'
    do i = 1, n
       write(unit, '(ES15.8)') stuff(i, 4)
    end do
 
    ! Close the file
    close(unit)
 
  end subroutine write_stuff_vtk
 
  subroutine func1(this, n, m, f0val, df0dx, fval, dfdx)
    ! ----------------------------------------------------------- !
    !  This subroutine calculates function values and gradients   !
    !  for "toy problem 3":                                       !
    !                                                             !
    !    minimize sum_(j = 1, .., n) xj/n                           !
    !  subject to sum_(j = 1, .., n) {(xj - Xj_GLL)^2} = 0          !
    ! ----------------------------------------------------------- !
    implicit none
    class(mma_comp_t), intent(inout) :: this
 
    integer, intent(in) :: n, m
    real(kind=rp), intent(inout) :: f0val
    real(kind=rp), dimension(n), intent(inout) :: df0dx
    real(kind=rp), dimension(m), intent(inout) :: fval
    real(kind=rp), dimension(m, n), intent(inout) :: dfdx
 
    real(kind=rp), dimension(n) :: x, coordx
    integer :: ierr, nglobal
    real(kind=rp) :: globalf0val
 
    call mpi_allreduce(n, nglobal, 1, &
         mpi_integer, mpi_sum, neko_comm, ierr)
 
    x = reshape(this%designx%x, [n])
    coordx = reshape(this%designx%dof%x, [n])
 
    f0val = sum(x) / nglobal
    df0dx = 1.0_rp / nglobal
    globalf0val = 0.0_rp
    call mpi_allreduce(f0val, globalf0val, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    f0val = globalf0val
    ! f0val = 0
    ! df0dx = 0
    fval(1) = sum((x-coordx)**2)
    globalf0val = 0.0_rp
    call mpi_allreduce(fval(1), globalf0val, 1, &
         mpi_real_precision, mpi_sum, neko_comm, ierr)
    fval(1) = globalf0val
 
    dfdx(1, :) = 2.0_rp * (x - coordx)
    fval(2) = - fval(1)
    dfdx(2, :) = - dfdx(1, :)
  end subroutine func1
 
end module mma_simcomp