FIXME : Add documentation
Pretty sure this barrier is not needed
call barrier
if (proc0) write(,) 'beginning class ',i,' with size ',nidxN_class(i)
Allocate matrix am. First dimension is basically theta along the entire
connected domain for each field. Second dimension is the local section
of the M_class(i)N_Class(i)(2ntgrid+1)nfield compound domain.
Clearly this will
allocate (am(nidxN_class(i), f_lo(i)%llim_proc:f_lo(i)%ulim_alloc))
Do we need to zero all 8 arrays on every loop? This can be more expensive than might think.
am = 0.0
call init_inverse_matrix (am, i)
Free memory
deallocate (am)
end do
DD> Comments
subroutine init_response_matrix
use mp, only: barrier, land_allreduce, proc0, nproc, iproc
use fields_arrays, only: phi, apar, bpar, phinew, aparnew, bparnew
use theta_grid, only: ntgrid
use kt_grids, only: naky, ntheta0
use dist_fn_arrays, only: g
use run_parameters, only: has_phi, has_apar, has_bpar
use gs2_layouts, only: init_fields_layouts, f_lo
use gs2_layouts, only: init_jfields_layouts
use file_utils, only: error_unit
use array_utils, only: zero_array
implicit none
integer :: ig, ifield, it, ik, i, m, n
complex, dimension(:,:), allocatable :: am
logical :: endpoint
logical :: response_was_read
nfield = 0
if (has_phi) nfield = nfield + 1
if (has_apar) nfield = nfield + 1
if (has_bpar) nfield = nfield + 1
nidx = (2*ntgrid+1)*nfield
call init_fields_layouts (nfield, nidx, naky, ntheta0, M_class, N_class, i_class, nproc, iproc)
call init_jfields_layouts (nfield, nidx, naky, ntheta0, i_class, nproc, iproc)
call finish_fields_layouts
response_was_read = .false.
!Either read the reponse
if(read_response) then
call read_response_from_file_imp(could_read = response_was_read)
! Ensure we reduce response_was_read across all processors so if _any_
! processors couldn't read the matrices we recalculate on all processors.
call land_allreduce(response_was_read)
if((.not. response_was_read) .and. proc0) write(error_unit(), &
'("Could not find response file so reverting to calculation.")')
!elseif(skip_initialisation) then
!do i = i_class, 1, -1
!!Pretty sure this barrier is not needed
!call barrier
!! if (proc0) write(*,*) 'beginning class ',i,' with size ',nidx*N_class(i)
!!Allocate matrix am. First dimension is basically theta along the entire
!!connected domain for each field. Second dimension is the local section
!!of the M_class(i)*N_Class(i)*(2*ntgrid+1)*nfield compound domain.
!!Clearly this will
!allocate (am(nidx*N_class(i), f_lo(i)%llim_proc:f_lo(i)%ulim_alloc))
!!Do we need to zero all 8 arrays on every loop? This can be more expensive than might think.
!am = 0.0
!call init_inverse_matrix (am, i)
!!Free memory
!deallocate (am)
!end do
end if
!or calculate it
if (.not. response_was_read) then
!
! keep storage cost down by doing one class at a time
! Note: could define a superclass (of all classes), a structure containing all am,
! then do this all at once. This would be faster, especially for large runs in a
! sheared domain, and could be triggered by local_field_solve
!
!<DD> Comments
!A class refers to a class of connected domain.
!These classes are defined by the extent of the connected domain, there can be
!many members of each class.
!There are i_class classes in total.
!N_class(ic) is a count of how many 2pi domains there are in members of class ic
!M_class(ic) is how many members of class ic there are.
!Sum N_class(ic)*M_class(ic) for ic=1,i_class is naky*ntheta0
!In comments cell refers to a 2pi domain whilst supercell is the connected domain,
!i.e. we have classes of supercells based on the number of cells they contain.
do i = i_class, 1, -1
!Pretty sure this barrier is not needed
call barrier
! if (proc0) write(*,*) 'beginning class ',i,' with size ',nidx*N_class(i)
!Allocate matrix am. First dimension is basically theta along the entire
!connected domain for each field. Second dimension is the local section
!of the M_class(i)*N_Class(i)*(2*ntgrid+1)*nfield compound domain.
!Clearly this will
allocate (am(nidx*N_class(i), f_lo(i)%llim_proc:f_lo(i)%ulim_alloc))
!Do we need to zero all 8 arrays on every loop? This can be more expensive than might think.
call zero_array(am)
call zero_array(g)
call zero_array(phi) ; call zero_array(phinew)
call zero_array(apar) ; call zero_array(aparnew)
call zero_array(bpar) ; call zero_array(bparnew)
!Loop over individual 2pi domains / cells
do n = 1, N_class(i)
!Loop over theta grid points in cell
!This is like a loop over nidx as we also handle all the fields in this loop
do ig = -ntgrid, ntgrid
!Are we at a connected boundary point on the lower side (i.e. left hand end of a
!tube/cell connected to the left)
endpoint = n > 1
endpoint = ig == -ntgrid .and. endpoint
!Start counting fields
ifield = 0
!Find response to phi
if (has_phi) then
ifield = ifield + 1
if (endpoint) then
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n-1)
it = f_lo(i)%it(m,n-1)
phinew(ntgrid,it,ik) = 1.0
end do
endif
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n)
it = f_lo(i)%it(m,n)
phinew(ig,it,ik) = 1.0
end do
if (.not. skip_initialisation) call init_response_row (ig, ifield, am, i, n)
call zero_array(phinew)
end if
!Find response to apar
if (has_apar) then
ifield = ifield + 1
if (endpoint) then
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n-1)
it = f_lo(i)%it(m,n-1)
aparnew(ntgrid,it,ik) = 1.0
end do
endif
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n)
it = f_lo(i)%it(m,n)
aparnew(ig,it,ik) = 1.0
end do
call init_response_row (ig, ifield, am, i, n)
call zero_array(aparnew)
end if
!Find response to bpar
if (has_bpar) then
ifield = ifield + 1
if (endpoint) then
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n-1)
it = f_lo(i)%it(m,n-1)
bparnew(ntgrid,it,ik) = 1.0
end do
endif
!Do all members of supercell together
do m = 1, M_class(i)
ik = f_lo(i)%ik(m,n)
it = f_lo(i)%it(m,n)
bparnew(ig,it,ik) = 1.0
end do
call init_response_row (ig, ifield, am, i, n)
call zero_array(bparnew)
end if
end do
end do
!Invert the matrix
call init_inverse_matrix (am, i)
!Free memory
deallocate (am)
end do
endif
if(dump_response) call dump_response_to_file_imp
end subroutine init_response_matrix