A routine to read the response matrix from file and populate the implicit response storage, note we also allocate the response storage objects DD>It may be anticipated that we need to fix the boundary points
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| logical, | intent(out) | :: | could_read | |||
| character(len=*), | intent(in), | optional | :: | suffix |
subroutine read_response_from_file_imp(could_read, suffix)
use fields_arrays, only: get_specific_response_file_name, time_read_response
use gs2_time, only: code_dt
use theta_grid, only: ntgrid
use kt_grids, only: naky, ntheta0, itright, get_leftmost_it
use gs2_layouts, only: f_lo, jf_lo, ij_idx, idx_local, idx, proc_id,idx_local
use mp, only: proc0, send, receive
use gs2_save, only: gs2_restore_response
use file_utils, only: error_unit
use job_manage, only: time_message
implicit none
logical, intent(out) :: could_read
character(len=*), optional, intent(in) :: suffix !If passed then use as part of file suffix
complex, dimension(:,:), allocatable :: tmp_arr, tmp_arr_full
complex, dimension(:), allocatable :: tmp_vec_full, tmp_vec
integer :: ic, im, ik, it, itmin, supercell_length, supercell_length_bound, in, ifld, ig, is_tmp
integer :: jflo, dc, nn, in_tmp, icount, it_tmp, nl, nr, ifld_tmp, ext_dom_length, ig_tmp, cur_idx
integer :: jflo_dup, dc_dup
integer, dimension(:,:), allocatable :: it_to_is, leftmost_it
integer, dimension(:), allocatable :: tmp_ints
logical :: is_local, is_local_dup
character(len = 256) :: file_name
logical :: file_found
real :: file_time_step
call time_message(.false.,time_read_response,' Field Read')
could_read = .true.
!First allocate the matrix storage
call alloc_response_objects
!Make a lookup array to convert itmin (the leftmost it in a connected domain)
!to the supercell index "is" used in the local fields. This will be used to
!ensure equivalent files can be given the same name.
allocate(it_to_is(ntheta0,naky),leftmost_it(ntheta0,naky),tmp_ints(ntheta0))
it_to_is=0
!//Note the following code is mostly borrowed from fm_init in the local fields
!First find all the leftmost it
do ik=1,naky
do it=1,ntheta0
leftmost_it(it,ik)=get_leftmost_it(it,ik)
enddo
enddo
!Now find supercell ids for each ky at a time
do ik=1,naky
tmp_ints=leftmost_it(:,ik)
it_tmp=0
is_tmp=0
do while(sum(tmp_ints)/=-1*ntheta0)
it_tmp=it_tmp+1
cur_idx=tmp_ints(it_tmp)
!If we've seen this domain skip
if(cur_idx==-1)cycle
!Increment counter
is_tmp=is_tmp+1
!Here we store the value
it_to_is(it_tmp,ik)=is_tmp
!Now we set all other connected locations to -1
!and store the appropriate is value
do it=1,ntheta0
if(tmp_ints(it)==cur_idx) then
tmp_ints(it)=-1
it_to_is(it,ik)=is_tmp
endif
enddo
enddo
enddo
!Cleanup
deallocate(tmp_ints)
!/End of borrowed code
!Notation recap:
! A class refers to all independent domains with the same length
! i_class is how many classes we have
! N_class(i_class) is how many 2Pi domains are in each member of i_class
! M_class(i_class) is how many independent domains are in i_class
allocate(tmp_vec(nfield*(2*ntgrid+1)))
allocate(tmp_arr(1+(2*ntgrid),nfield))
!Loop over classes (supercell length)
do ic=1,i_class
!Extended domain length
ext_dom_length=1+(2*ntgrid)*N_class(ic)
!Work out how long the supercell is
supercell_length=nfield*ext_dom_length !Without boundary points
supercell_length_bound=(1+2*ntgrid)*nfield*N_class(ic) !With boundary points
!Make storage
allocate(tmp_arr_full(supercell_length,supercell_length))
allocate(tmp_vec_full(supercell_length))
!Now loop over all members of this class
do im=1,M_class(ic)
tmp_arr_full=0.
tmp_vec_full=0.
!Now we are thinking about a single supercell
!we can get certain properties before looping
!over the individual elements
!Get the ik index
ik=f_lo(ic)%ik(im,1)
!Get the leftmost it index (named itmin to match local field routines)
!This is currently used to identify the supercell like "is" is used in
!the local field routines. It would be nice to also use "is" here (or
!"itmin" there).
itmin=leftmost_it(f_lo(ic)%it(im,1),ik)
!Now make file name
file_name = adjustl(get_specific_response_file_name(ik, it_to_is(itmin, ik), code_dt, suffix))
! First check if the expected file exists, if not then exit
inquire( file = trim(file_name), exist = file_found)
if (.not. file_found) then
if(proc0) write(error_unit(), &
'("Could not find response file ",A,", reverting to calculation.")') trim(file_name)
could_read = .false.
return
end if
if(proc0)then
call gs2_restore_response(tmp_arr_full, file_name, file_time_step)
! Should perhaps double check that file_time_step == code_dt?
endif
!Initialise counter
icount=1
!Loop over the fields
do ifld=1,nfield
!Reinitialise "it" back to the start of the supercell chain
it=itmin
!Loop over the different it (2Pi domains)
do in=1,N_class(ic)
!Loop over theta
do ig=-ntgrid,ntgrid
!Skip the duplicate boundary points -- This is no good here. !<DD>
! if((ig.eq.ntgrid).and.(in.ne.N_class(ic))) cycle
!Convert to jf_lo index
jflo=ij_idx(jf_lo,ig,ifld,ik,it)
!See if it's local
is_local=idx_local(jf_lo,jflo)
!If it's not local then we have nothing to do
!unless we're the proc who writes (proc0).
if(.not.(is_local.or.proc0)) cycle
!Get row
if(proc0)then
tmp_vec_full=tmp_arr_full(:,icount)
!Increment counter
if(.not.(ig==ntgrid.and.in/=N_Class(ic))) icount=icount+1
endif
!Now unpack tmp_vec_full and do communications if needed
if(is_local)then
!No comms needed if local
if(.not.proc0) call receive(tmp_vec_full,0)
!Get dcell index
dc=jf_lo%dj(ic,jflo)
!Now we pack the tmp_vec in the correct order
!We must fill in the boundary points
!We need to pick the value of "n" in the right order
it_tmp=itmin
do in_tmp=1,N_class(ic)
tmp_arr=0
tmp_vec=0
!Now we need to work out where to put things in tmp_vec_full
!In doing this we try to match the local fields data layout
!to aid comparisons
do ifld_tmp=1,nfield
do ig_tmp=1,2*ntgrid+1
!Skip boundary points
if((ig_tmp==(2*ntgrid+1)).and.(in_tmp/=N_class(ic))) cycle
!Get index
cur_idx=ig_tmp+(2*ntgrid)*(in_tmp-1)+(ifld_tmp-1)*ext_dom_length
!Store data
tmp_arr(ig_tmp,ifld_tmp)=tmp_vec_full(cur_idx)
enddo
enddo
!<DD>It may be anticipated that we need to fix the boundary points
!here but we don't actually need to do anything.
!Because we sum over the entire supercell in getfield we only want
!the repeated boundary point to be included once.
!We still need to calculate the field at the repeated point but the
!fix for that is handled at the bottom of the routine
!In other words we don't need something of the form:
! !Fix boundary points
! if(in_tmp.ne.N_class(ic))then
! do ifld_tmp=1,nfield
! cur_idx=1+(2*ntgrid)*(in_tmp)+(ifld_tmp-1)*ext_dom_length
! tmp_arr(2*ntgrid+1,ifld_tmp)=tmp_vec_full(cur_idx)
! enddo
! endif
!Store in correct order
do ifld_tmp=1,nfield
nl=1+(ifld_tmp-1)*(2*ntgrid+1)
nr=nl+2*ntgrid
tmp_vec(nl:nr)=tmp_arr(:,ifld_tmp)
enddo
!Pick the correct n
do nn=1,N_class(ic)
if(f_lo(ic)%it(im,nn)==it_tmp) exit
enddo
!Now we can get supercell range (including boundaries)
nl=1+nidx*(nn-1)
nr=nl+nidx-1
!Store section
aminv(ic)%dcell(dc)%supercell(nl:nr)=tmp_vec
!If we're not at the last cell then increment it_tmp.
!This check is needed to avoid trying to access itright
!in non-box simulations where this array is not allocated.
if(in_tmp < N_class(ic)) then
it_tmp=itright(it_tmp, ik)
end if
enddo
else
!Only proc0 should get here but test anyway
if(proc0) call send(tmp_vec_full,proc_id(jf_lo,jflo))
endif
enddo
!If we're not at the last cell then increment it.
!This check is needed to avoid trying to access itright
!in non-box simulations where this array is not allocated.
if(in < N_class(ic)) then
it=itright(it, ik)
end if
enddo
enddo
!Now we need to fill in the repeated boundary points
!If there are no boundary points then advance
if(N_class(ic)==1) cycle
it=itmin
do in=1,N_class(ic)-1
do ifld=1,nfield
!First get the index of the point we want to fill
jflo=ij_idx(jf_lo,ntgrid,ifld,ik,it)
!Now we get the index of the point which has this data
jflo_dup=ij_idx(jf_lo,-ntgrid,ifld,ik,itright(it, ik))
!Now get locality
is_local=idx_local(jf_lo,jflo)
is_local_dup=idx_local(jf_lo,jflo_dup)
!Get dcell values
if(is_local) dc=jf_lo%dj(ic,jflo)
if(is_local_dup) dc_dup=jf_lo%dj(ic,jflo_dup)
!Now copy/communicate
if(is_local)then
if(is_local_dup)then
aminv(ic)%dcell(dc)%supercell=aminv(ic)%dcell(dc_dup)%supercell
else
call receive(aminv(ic)%dcell(dc)%supercell,proc_id(jf_lo,jflo_dup))
endif
elseif(is_local_dup)then
call send(aminv(ic)%dcell(dc_dup)%supercell,proc_id(jf_lo,jflo))
endif
enddo
!Increment it -- note we don't need to guard against access to itright
!here as for non-box runs we won't enter the outer in=1,N_class(ic)-1
!loop as N_class(ic) == 1
it=itright(it,ik)
enddo
end do
!Free
deallocate(tmp_arr_full,tmp_vec_full)
end do
!Tidy
deallocate(tmp_vec,tmp_arr,leftmost_it,it_to_is)
call time_message(.false.,time_read_response,' Field Read')
end subroutine read_response_from_file_imp