! simulation module for QuickPIC Open Source 1.0 ! update: 01/09/2018 module simulation_class use parallel_class use parallel_pipe_class use perrors_class use spect3d_class use spect2d_class use field2d_class use fdist2d_class use field3d_class use beam3d_class use species2d_class use fdist3d_class use hdf5io_class use input_class use mpi implicit none private public :: simulation type fdist2d_wrap class(fdist2d), allocatable :: p end type fdist2d_wrap ! type fdist3d_wrap class(fdist3d), allocatable :: p end type fdist3d_wrap ! type sim_fields private class(parallel_pipe),pointer :: p => null() class(perrors),pointer :: err => null() class(spect3d), pointer :: sp3 => null() class(spect2d), pointer :: sp2 => null() type(field2d), allocatable :: qb, qe, psit, psi, div_vpot, reg type(field2d), allocatable :: fxy, bxyz, cu, dcu, amu, epw, epwb type(field3d), allocatable :: bexyz, bbxyz type(field3d), allocatable :: psi3d,cu3d contains generic :: new => init_sim_fields generic :: del => end_sim_fields procedure, private :: init_sim_fields, end_sim_fields end type sim_fields ! type sim_beams private class(parallel_pipe),pointer :: p => null() class(perrors),pointer :: err => null() class(spect3d), pointer :: sp3 => null() class(spect2d), pointer :: sp2 => null() type(beam3d), dimension(:), allocatable :: beam type(fdist3d_wrap), dimension(:), allocatable :: pf contains generic :: new => init_sim_beams generic :: del => end_sim_beams procedure, private :: init_sim_beams, end_sim_beams end type sim_beams ! type sim_species private class(parallel_pipe),pointer :: p => null() class(perrors),pointer :: err => null() class(spect3d), pointer :: sp3 => null() class(spect2d), pointer :: sp2 => null() type(fdist2d_wrap), dimension(:), allocatable :: pf type(species2d), dimension(:), allocatable :: spe contains generic :: new => init_sim_species generic :: del => end_sim_species procedure, private :: init_sim_species, end_sim_species end type sim_species ! type sim_diag private class(parallel_pipe),pointer :: p => null() class(perrors),pointer :: err => null() class(spect3d), pointer :: sp3 => null() class(spect2d), pointer :: sp2 => null() type(hdf5file) :: file class(*), pointer :: obj => null() integer, allocatable :: slice, slice_pos, psample, dim integer :: df end type sim_diag ! type simulation private ! type(input) :: sim type(input_json), pointer :: in => null() class(parallel_pipe),pointer :: p => null() class(perrors),pointer :: err => null() class(spect3d), pointer :: sp3 => null() class(spect2d), pointer :: sp2 => null() type(sim_fields) :: fields type(sim_beams) :: beams type(sim_species) :: species type(sim_diag), dimension(:), allocatable :: diag integer :: iter, nstep3d, nstep2d, start3d,nbeams,nspecies, tstep integer, dimension(8) :: tag integer, dimension(:), allocatable :: tag_spe, id_spe, id integer, dimension(:), allocatable :: tag_beam, id_beam integer, dimension(:,:), allocatable :: id_bq, tag_bq real :: dex, dxi, dex2, dt contains generic :: new => init_simulation generic :: del => end_simulation generic :: go => go_simulation procedure, private :: init_simulation, end_simulation procedure, private :: init_diag, diag_simulation procedure, private :: go_simulation end type simulation ! character(len=20), save :: class = 'simulation: ' character(len=128), save :: erstr contains ! subroutine init_simulation(this) implicit none class(simulation), intent(inout) :: this ! local data character(len=18), save :: sname = 'init_simulation:' real :: min, max, n0, dx, dy, dz, dt integer :: indx, indy, indz logical :: read_rst allocate(this%in) call this%in%new() this%err => this%in%err this%p => this%in%pp this%sp3 => this%in%sp this%sp2 => this%in%sp call this%err%werrfl2(class//sname//' started') call this%in%get('simulation.n0',n0) call this%in%get('simulation.indx',indx) call this%in%get('simulation.indy',indy) call this%in%get('simulation.indz',indz) call this%in%get('simulation.box.x(1)',min) call this%in%get('simulation.box.x(2)',max) dx=(max-min)/real(2**indx) call this%in%get('simulation.box.y(1)',min) call this%in%get('simulation.box.y(2)',max) dy=(max-min)/real(2**indy) call this%in%get('simulation.box.z(1)',min) call this%in%get('simulation.box.z(2)',max) dz=(max-min)/real(2**indz) this%dex = dx this%dxi = dz this%dex2 = dx * dx this%nstep2d = 2**this%sp3%getindz()/this%p%getnstage() call this%in%get('simulation.time',min) call this%in%get('simulation.dt',max) this%nstep3d = min/max call this%in%get('simulation.read_restart',read_rst) call this%in%get('simulation.dt',dt) this%dt = dt if (read_rst) then call this%in%get('simulation.restart_timestep',this%start3d) this%start3d = this%start3d + 1 else this%start3d = 1 end if call this%in%get('simulation.iter',this%iter) call this%in%get('simulation.nbeams',this%nbeams) call this%in%get('simulation.nspecies',this%nspecies) call this%fields%new(this%in) call this%beams%new(this%in,this%fields) call this%species%new(this%in,this%fields,(this%start3d-1)*dt) call this%init_diag() allocate(this%tag_spe(this%nspecies),this%tag_beam(this%nbeams)) allocate(this%id_spe(this%nspecies),this%id_beam(this%nbeams)) allocate(this%id_bq(this%nbeams,3),this%tag_bq(this%nbeams,2)) allocate(this%id(9+size(this%diag))) this%id(:) = MPI_REQUEST_NULL this%id_spe(:) = MPI_REQUEST_NULL this%id_beam(:) = MPI_REQUEST_NULL this%id_bq(:,:) = MPI_REQUEST_NULL call this%err%werrfl2(class//sname//' ended') end subroutine init_simulation ! subroutine end_simulation(this) implicit none class(simulation), intent(inout) :: this ! local data character(len=18), save :: sname = 'end_simulation:' integer :: ierr call this%err%werrfl2(class//sname//' started') ! call this%fields%del() ! call this%beams%del() ! call this%species%del() call MPI_FINALIZE(ierr) call this%err%werrfl2(class//sname//' ended') end subroutine end_simulation ! subroutine init_sim_fields(this,input) implicit none class(sim_fields), intent(inout) :: this type(input_json), pointer, intent(inout) :: input ! local data character(len=18), save :: sname = 'init_sim_fields:' character(len=18), save :: class = 'sim_fields:' character(len=20) :: s1, s2, s3 character(len=:), allocatable :: ff integer :: i,n,ndump,j,k,l,m this%err => input%err this%p => input%pp this%sp3 => input%sp this%sp2 => input%sp call this%err%werrfl2(class//sname//' started') allocate(this%qb, this%qe, this%psit, this%psi) allocate(this%div_vpot, this%reg, this%fxy, this%bxyz, this%cu) allocate(this%dcu, this%amu, this%epw, this%epwb) allocate(this%bexyz, this%bbxyz) call this%bexyz%new(this%p,this%err,this%sp3,dim=3) call this%bbxyz%new(this%p,this%err,this%sp3,dim=3) call this%qb%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.,gcells=1) call this%qe%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.) call this%psit%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.) call this%div_vpot%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.) call this%psi%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.) call this%reg%new(this%p,this%err,this%sp2,dim=1,fftflag=.true.) call this%fxy%new(this%p,this%err,this%sp2,dim=2,fftflag=.true.) call this%cu%new(this%p,this%err,this%sp2,dim=3,fftflag=.true.,state=1) call this%dcu%new(this%p,this%err,this%sp2,dim=2,fftflag=.true.) call this%amu%new(this%p,this%err,this%sp2,dim=3,fftflag=.true.) call this%epw%new(this%p,this%err,this%sp2,dim=2,fftflag=.true.,state=1) call this%epwb%new(this%p,this%err,this%sp2,dim=2,fftflag=.true.) call this%bxyz%new(this%p,this%err,this%sp2,dim=3,fftflag=.true.) call input%get('simulation.nspecies',n) loop1: do i = 1, n write (s1, '(I4.4)') i call input%info('species('//trim(s1)//').diag',n_children=m) do j = 1, m write (s2, '(I4.4)') j call input%get('species('//trim(s1)//').diag'//'('//trim(s2)//').ndump',ndump) if (ndump>0) then call input%info('species('//trim(s1)//').diag'//'('//trim(s2)//').name',n_children=l) do k = 1, l write (s3, '(I4.4)') k if(allocated(ff)) deallocate(ff) call input%get('species('//trim(s1)//').diag'//'('//trim(s2)//').name'& &//'('//trim(s3)//')',ff) if (ff == 'jx' .or. ff == 'jy' .or. ff == 'jz') then allocate(this%cu3d) call this%cu3d%new(this%p,this%err,this%sp3,dim=1) exit loop1 end if end do end if end do end do loop1 call input%info('field.diag',n_children=n) loop2: do i = 1, n write (s1,'(I4.4)') i call input%get('field.diag('//trim(s1)//').ndump',ndump) if (ndump > 0) then call input%info('field.diag('//trim(s1)//').name',n_children=m) do j = 1, m write (s2,'(I4.4)') j if(allocated(ff)) deallocate(ff) call input%get('field.diag('//trim(s1)//').name('//trim(s2)//')',ff) if (ff == 'psi') then allocate(this%psi3d) call this%psi3d%new(this%p,this%err,this%sp3,dim=1) exit loop2 end if end do end if end do loop2 call this%err%werrfl2(class//sname//' ended') end subroutine init_sim_fields ! subroutine end_sim_fields(this) implicit none class(sim_fields), intent(inout) :: this ! local data character(len=18), save :: sname = 'end_sim_fields:' character(len=18), save :: class = 'sim_fields:' integer :: i, n call this%err%werrfl2(class//sname//' started') call this%bexyz%del() call this%bbxyz%del() call this%qb%del() call this%qe%del() call this%psit%del() call this%div_vpot%del() call this%psi%del() call this%reg%del() call this%fxy%del() call this%cu%del() call this%dcu%del() call this%amu%del() call this%epw%del() call this%epwb%del() call this%bxyz%del() if(allocated(this%psi3d)) call this%psi3d%del() if(allocated(this%cu3d)) call this%cu3d%del() call this%err%werrfl2(class//sname//' ended') end subroutine end_sim_fields ! subroutine init_sim_beams(this,input,fields) implicit none class(sim_beams), intent(inout) :: this type(input_json), pointer, intent(inout) :: input class(sim_fields), intent(inout) :: fields ! local data character(len=18), save :: class = 'sim_beams:' character(len=18), save :: sname = 'init_sim_beams:' integer :: i,n integer :: npf,npx,npy,npz,npp real, dimension(3,100) :: arg logical :: quiet real :: gamma real :: qm, qbm, dt logical :: read_rst integer :: rst_timestep, ierr type(hdf5file) :: file_rst character(len=20) :: sn, sid, stime,s1 this%err => input%err this%p => input%pp this%sp3 => input%sp this%sp2 => input%sp call this%err%werrfl2(class//sname//' started') call input%get('simulation.nbeams',n) allocate(this%beam(n),this%pf(n)) do i = 1, n arg(:,:) = 0.0 write (sn,'(I3.3)') i s1 = 'beam('//trim(sn)//')' call input%get(trim(s1)//'.profile',npf) select case (npf) case (0) allocate(fdist3d_000::this%pf(i)%p) call this%pf(i)%p%new(input,i) ! Add new distributions under this line case default write (erstr,*) 'Invalid species profile number:', npf call this%err%equit(class//sname//erstr) end select call input%get(trim(s1)//'.q',qm) call input%get(trim(s1)//'.m',qbm) qbm = qm/qbm call input%get('simulation.dt',dt) call this%beam(i)%new(this%p,this%err,this%sp3,this%pf(i)%p,qbm=qbm,& &dt=dt,ci=1.0,xdim=7) call input%get('simulation.read_restart',read_rst) if (read_rst) then call input%get('simulation.restart_timestep',rst_timestep) write (sn,'(I4.4)') i write (sid,'(I10.10)') this%p%getidproc() write (stime,'(I8.8)') rst_timestep call file_rst%new(& &filename = './RST/Beam-'//trim(sn)//'/',& &dataname = 'RST-beam'//trim(sn)//'-'//trim(sid),& &n = rst_timestep) call this%beam(i)%rrst(file_rst) end if end do call MPI_BARRIER(this%p%getlworld(),ierr) call this%err%werrfl2(class//sname//' ended') end subroutine init_sim_beams ! subroutine end_sim_beams(this) implicit none class(sim_beams), intent(inout) :: this ! local data character(len=18), save :: sname = 'end_sim_beams:' character(len=18), save :: class = 'sim_beams:' integer :: i, n call this%err%werrfl2(class//sname//' started') n = size(this%beam) do i = 1, n call this%beam(i)%del() end do call this%err%werrfl2(class//sname//' ended') end subroutine end_sim_beams ! subroutine init_sim_species(this,input,fields,s) implicit none class(sim_species), intent(inout) :: this type(input_json), pointer, intent(inout) :: input class(sim_fields), intent(inout) :: fields real, intent(in) :: s ! local data character(len=18), save :: class = 'sim_species:' character(len=18), save :: sname = 'init_sim_species:' integer :: i,n,ndump integer :: npf real, dimension(3,100) :: arg character(len=20) :: sn,s1 integer :: indz, xppc, yppc real :: min, max, cwp, n0 real :: qm, qbm, dz this%err => input%err this%p => input%pp this%sp3 => input%sp this%sp2 => input%sp call this%err%werrfl2(class//sname//' started') call input%get('simulation.n0',n0) call input%get('simulation.indz',indz) cwp=5.32150254*1e9/sqrt(n0) call input%get('simulation.box.z(1)',min) call input%get('simulation.box.z(2)',max) dz=(max-min)/real(2**indz) call input%get('simulation.nspecies',n) allocate(this%spe(n),this%pf(n)) do i = 1, n write (sn,'(I3.3)') i s1 = 'species('//trim(sn)//')' call input%get(trim(s1)//'.profile',npf) select case (npf) case (0) allocate(fdist2d_000::this%pf(i)%p) call this%pf(i)%p%new(input,i) ! Add new distributions under this line case default write (erstr,*) 'Invalid species profile number:', npf call this%err%equit(class//sname//erstr) end select call input%get(trim(s1)//'.q',qm) call input%get(trim(s1)//'.m',qbm) qbm = qm/qbm call this%spe(i)%new(this%p,this%err,this%sp3,this%pf(i)%p,& &qbm=qbm,dt=dz,ci=1.0,xdim=8,s=s) end do call this%err%werrfl2(class//sname//' ended') end subroutine init_sim_species ! subroutine end_sim_species(this) implicit none class(sim_species), intent(inout) :: this ! local data character(len=18), save :: sname = 'end_sim_species:' character(len=18), save :: class = 'sim_species:' integer :: i, n call this%err%werrfl2(class//sname//' started') n = size(this%spe) do i = 1, n call this%spe(i)%del() end do call this%err%werrfl2(class//sname//' ended') end subroutine end_sim_species ! subroutine go_simulation(this) implicit none class(simulation), intent(inout) :: this ! local data character(len=18), save :: sname = 'go_simulation:' integer, dimension(10) :: istat integer :: i,j,l,m,n,ierr call this%err%werrfl2(class//sname//' started') do i = this%start3d, this%nstep3d this%tstep = i write (erstr,*) '3D step:', i call this%err%werrfl0(erstr) do m = 1, this%nbeams this%tag_bq(m,1) = ntag() this%tag_bq(m,2) = ntag() call this%beams%beam(m)%qdp(this%id_bq(m,1),this%id_bq(m,2),& &this%id_bq(m,3),this%tag_bq(m,1),this%tag_bq(m,2)) end do do l = 1, this%nspecies this%tag_spe(l) = ntag() call this%species%spe(l)%precv(this%tag_spe(l)) end do this%tag(3) = ntag() call this%fields%cu%precv(this%tag(3)) this%tag(4) = ntag() call this%fields%epw%precv(this%tag(4)) this%tag(5) = ntag() call this%fields%fxy%precv(this%tag(5)) this%tag(6) = ntag() call this%fields%psit%precv(this%tag(6)) this%tag(7) = ntag() call this%fields%bxyz%precv(this%tag(7)) call this%fields%fxy%cb(this%fields%bexyz,1,(/1,2/),(/1,2/)) call this%fields%psit%cb(this%fields%bexyz,1,(/1/),(/3/)) call this%fields%bxyz%cb(this%fields%bbxyz,1,(/1,2,3/),(/1,2,3/)) do j = 1, this%nstep2d write (erstr,*) '2D step:', j call this%err%werrfl0(erstr) if (j == this%nstep2d) then do m = 1, this%nbeams call MPI_WAIT(this%id_bq(m,2),istat,ierr) end do endif call this%fields%qb%as(0.0) do m = 1, this%nbeams call this%beams%beam(m)%qdp(this%fields%qb,j+1) end do call this%fields%qb%fftrk(1) call this%fields%qb%elf(this%fields%epwb) call this%fields%qe%as(0.0) do l = 1, this%nspecies call this%species%spe(l)%qdp(this%fields%qe) end do call this%fields%qe%fftrk(1) call this%fields%qe%pot(this%fields%psi) call this%fields%psi%grad(this%fields%fxy) call this%fields%fxy%fftkr(1) call this%fields%psi%fftkr(1) do l = 1, this%nspecies call this%species%spe(l)%extpsi(this%fields%psi,this%dex) end do call this%fields%psi%mult(this%fields%psi,this%dex*this%dex) if (allocated(this%fields%psi3d)) call this%fields%psi%cb(this%fields%psi3d,j+1,(/1/),(/1/)) call this%fields%cu%div(this%fields%div_vpot) call this%fields%div_vpot%pot(this%fields%psit) call this%fields%psit%mult(this%fields%psit,-this%dex) call this%fields%psit%fftkr(1) call this%fields%cu%bf(this%fields%bxyz) do l = 1, this%iter call this%fields%bxyz%mult(this%fields%epw,(/1,2/),(/2,1/),(/-this%dex,this%dex/)) call this%fields%bxyz%mult(this%fields%bxyz,(/3/),(/3/),(/this%dex/)) call this%fields%bxyz%sub(this%fields%bxyz,this%fields%epwb,(/1/),(/1/),(/2/)) call this%fields%bxyz%add(this%fields%bxyz,this%fields%epwb,(/2/),(/2/),(/1/)) call this%fields%bxyz%fftkr(2) call this%fields%cu%as(0.0) call this%fields%dcu%as(0.0) call this%fields%amu%as(0.0) do m = 1, this%nspecies call this%species%spe(m)%amjdp(this%fields%fxy,this%fields%bxyz,this%fields%psit,& &this%fields%cu,this%fields%amu,this%fields%dcu,this%dex) end do call this%fields%cu%mult(this%fields%cu,this%dex) call this%fields%amu%mult(this%fields%amu,this%dex) call this%fields%dcu%mult(this%fields%dcu,this%dex) if (l == this%iter) then do m = 1, this%nspecies call this%species%spe(m)%cbq(j+1) end do if (allocated(this%fields%cu3d)) then call this%fields%cu%cb(this%fields%cu3d,j+1,(/1,2,3/),(/1,2,3/)) end if endif call this%fields%cu%fftrk(1) call this%fields%dcu%fftrk(1) call this%fields%amu%fftrk(3) call this%fields%cu%bfqp(this%fields%dcu,this%fields%amu,this%fields%epw,this%dex2,& &this%dex) call this%fields%cu%div(this%fields%div_vpot) call this%fields%div_vpot%pot(this%fields%psit) call this%fields%psit%mult(this%fields%psit,-this%dex) call this%fields%psit%fftkr(1) call this%fields%cu%bf(this%fields%bxyz) enddo call this%fields%bxyz%mult(this%fields%epw,(/1,2/),(/2,1/),(/-this%dex,this%dex/)) call this%fields%bxyz%mult(this%fields%bxyz,(/3/),(/3/),(/this%dex/)) call this%fields%bxyz%sub(this%fields%bxyz,this%fields%epwb,(/1/),(/1/),(/2/)) call this%fields%bxyz%add(this%fields%bxyz,this%fields%epwb,(/2/),(/2/),(/1/)) call this%fields%bxyz%fftkr(2) call this%fields%fxy%mult(this%fields%fxy,-1.0) call this%fields%fxy%add(this%fields%fxy,this%fields%bxyz,(/1/),(/1/),(/2/)) call this%fields%fxy%sub(this%fields%fxy,this%fields%bxyz,(/2/),(/2/),(/1/)) call this%fields%dcu%mult(this%fields%dcu,this%dxi) call this%fields%cu%sub(this%fields%cu,this%fields%dcu,(/1,2/),(/1,2/),(/1,2/)) do m = 1, this%nspecies call this%species%spe(m)%push(this%fields%fxy,this%fields%bxyz,this%fields%psit,& &this%dex) end do call this%fields%fxy%mult(this%fields%fxy,this%dex) call this%fields%fxy%cb(this%fields%bexyz,j+1,(/1,2/),(/1,2/)) call this%fields%psit%cb(this%fields%bexyz,j+1,(/1/),(/3/)) call this%fields%bxyz%mult(this%fields%bxyz,(/1,2/),(/1,2/),(/this%dex,this%dex/)) call this%fields%bxyz%cb(this%fields%bbxyz,j+1,(/1,2,3/),(/1,2,3/)) enddo do m = 1, this%nspecies call this%species%spe(m)%psend(this%tag_spe(m),this%id_spe(m)) end do call MPI_WAIT(this%id(5),istat,ierr) call this%fields%cu%psend(this%tag(3),this%id(5)) call MPI_WAIT(this%id(6),istat,ierr) call this%fields%epw%psend(this%tag(4),this%id(6)) call MPI_WAIT(this%id(7),istat,ierr) call this%fields%fxy%psend(this%tag(5),this%id(7)) call MPI_WAIT(this%id(8),istat,ierr) call this%fields%psit%psend(this%tag(6),this%id(8)) call MPI_WAIT(this%id(9),istat,ierr) call this%fields%bxyz%psend(this%tag(7),this%id(9)) do m = 1, this%nbeams this%tag_beam(m) = ntag() call MPI_WAIT(this%id_beam(m),istat,ierr) call this%beams%beam(m)%push(this%fields%bexyz,this%fields%bbxyz,this%dex,this%dxi,& &this%tag_beam(m),this%tag_beam(m),this%id_beam(m)) end do call this%diag_simulation() do m = 1, this%nspecies call MPI_WAIT(this%id_spe(m),istat,ierr) call this%species%spe(m)%renew(i*this%dt) end do end do call this%err%werrfl2(class//sname//' ended') end subroutine go_simulation ! subroutine init_diag(this) implicit none class(simulation), intent(inout), target :: this ! local data character(len=18), save :: sname = 'init_diag:' integer :: n_diag = 0, ndump, slice, slice_pos, sample integer :: n, m, l, i, j, k, ii character(len=20) :: s1, s2, s3, s4, sn1, sn2, sn3, sn4 character(len=:), allocatable :: ss,sl real :: min, max, n0, dt real :: alx1, aly1, alz1, alx2, aly2, alz2 logical :: rst integer :: ierr, indx, indy, indz, dim class (*), pointer :: obj => null() call this%err%werrfl2(class//sname//' started') call this%in%get('simulation.n0',n0) call this%in%get('simulation.indx',indx) call this%in%get('simulation.indy',indy) call this%in%get('simulation.indz',indz) call this%in%get('simulation.box.x(1)',min) call this%in%get('simulation.box.x(2)',max) alx1 = min alx2 = max call this%in%get('simulation.box.y(1)',min) call this%in%get('simulation.box.y(2)',max) aly1 = min aly2 = max call this%in%get('simulation.box.z(1)',min) call this%in%get('simulation.box.z(2)',max) alz1 = min alz2 = max call this%in%get('simulation.dt',dt) do i = 1, this%nbeams write (s1, '(I4.4)') i call this%in%info('beam('//trim(s1)//').diag',n_children=m) do j = 1, m write (s2, '(I4.4)') j call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').ndump',ndump) if (ndump>0) then call this%in%info('beam('//trim(s1)//').diag'//'('//trim(s2)//').name',n_children=l) if(this%in%found('beam('//trim(s1)//').diag'//'('//trim(s2)//').slice')) then call this%in%info('beam('//trim(s1)//').diag'//'('//trim(s2)//').slice',n_children=n) n_diag = n_diag + l*n else n_diag = n_diag + l end if end if end do end do do i = 1, this%nspecies write (s1, '(I4.4)') i call this%in%info('species('//trim(s1)//').diag',n_children=m) do j = 1, m write (s2, '(I4.4)') j call this%in%get('species('//trim(s1)//').diag'//'('//trim(s2)//').ndump',ndump) if (ndump>0) then call this%in%info('species('//trim(s1)//').diag'//'('//trim(s2)//').name',n_children=l) if(this%in%found('species('//trim(s1)//').diag'//'('//trim(s2)//').slice')) then call this%in%info('species('//trim(s1)//').diag'//'('//trim(s2)//').slice',n_children=n) n_diag = n_diag + l*n else n_diag = n_diag + l end if end if end do end do call this%in%info('field.diag',n_children=n) do i = 1, n write (s1, '(I4.4)') i call this%in%get('field.diag('//trim(s1)//').ndump',ndump) if (ndump>0) then call this%in%info('field.diag('//trim(s1)//').name',n_children=l) if(this%in%found('field.diag('//trim(s1)//').slice')) then call this%in%info('field.diag('//trim(s1)//').slice',n_children=m) n_diag = n_diag + l*m else n_diag = n_diag + l end if end if end do call this%in%get('simulation.dump_restart',rst) if (rst) n_diag = n_diag + this%nbeams allocate(this%diag(n_diag)) n_diag = 0 do i = 1, this%nbeams write (s1, '(I4.4)') i call this%in%info('beam('//trim(s1)//').diag',n_children=m) do j = 1, m write (s2, '(I4.4)') j call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').ndump',ndump) if (ndump>0) then call this%in%info('beam('//trim(s1)//').diag'//'('//trim(s2)//').name',n_children=l) do k = 1, l write (s3, '(I4.4)') k if (allocated(ss)) deallocate(ss) call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').name'& &//'('//trim(s3)//')',ss) select case (trim(ss)) case ('charge') if (this%in%found('beam('//trim(s1)//').diag'//'('//trim(s2)//').slice')) then call this%in%info('beam('//trim(s1)//').diag'//'('//trim(s2)//').slice',n_children=n) do ii = 1, n n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => this%beams%beam(i) allocate(this%diag(n_diag)%slice,this%diag(n_diag)%slice_pos) allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = 1 if (allocated(sl)) deallocate(sl) write (s4, '(I4.4)') ii call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').& &slice('//trim(s4)//').(1)',sl) select case (sl) case ('yz') this%diag(n_diag)%slice = 1 call this%diag(n_diag)%file%new(& &axismin = (/aly1,alz1,alx1/),& &axismax = (/aly2,alz2,alx2/),& &axisname = (/'y ','\xi','x '/),& &axislabel = (/'y ','\xi','x '/)) case ('xz') this%diag(n_diag)%slice = 2 call this%diag(n_diag)%file%new(& &axismin = (/alx1,alz1,aly1/),& &axismax = (/alx2,alz2,aly2/),& &axisname = (/'x ','\xi','y '/),& &axislabel = (/'x ','\xi','y '/)) case ('xy') this%diag(n_diag)%slice = 3 call this%diag(n_diag)%file%new(& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &axisname = (/'x','y','z'/),& &axislabel = (/'x','y','z'/)) end select call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').& &slice('//trim(s4)//').(2)',this%diag(n_diag)%slice_pos) if (this%p%getidproc() == 0) then call system('mkdir -p ./Beam'//trim(s1)//'/Charge_slice_'//trim(s4)//'/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &rank = 2,& &filename = './Beam'//trim(s1)//'/Charge_slice_'//trim(s4)//'/',& &dataname = 'charge_slice_'//sl,& &units = 'n_0',& &label = 'Charge Density') end do else n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => this%beams%beam(i) allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = 1 if (this%p%getidproc() == 0) then call system('mkdir -p ./Beam'//trim(s1)//'/Charge/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisname = (/'x ','y ','\xi'/),& &axislabel = (/'x ','y ','\xi'/),& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &rank = 3,& &filename = './Beam'//trim(s1)//'/Charge/',& &dataname = 'charge',& &units = 'n_0',& &label = 'Charge Density') end if case ('raw') n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => this%beams%beam(i) allocate(this%diag(n_diag)%psample) call this%in%get('beam('//trim(s1)//').diag'//'('//trim(s2)//').& &sample',this%diag(n_diag)%psample) if (this%p%getidproc() == 0) then call system('mkdir -p ./Beam'//trim(s1)//'/Raw/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &ty = 'particles',& &filename = './Beam'//trim(s1)//'/Raw/',& &dataname = 'raw',& &units = '',& &label = 'Beam Raw') end select end do end if end do end do do i = 1, this%nspecies write (s1, '(I4.4)') i call this%in%info('species('//trim(s1)//').diag',n_children=m) do j = 1, m write (s2, '(I4.4)') j call this%in%get('species('//trim(s1)//').diag'//'('//trim(s2)//').ndump',ndump) if (ndump>0) then call this%in%info('species('//trim(s1)//').diag'//'('//trim(s2)//').name',n_children=l) do k = 1, l write (s3, '(I4.4)') k if (allocated(ss)) deallocate(ss) call this%in%get('species('//trim(s1)//').diag'//'('//trim(s2)//').name'& &//'('//trim(s3)//')',ss) select case (trim(ss)) case ('charge') sn1 = 'Charge' sn2 = 'charge' sn3 = 'n_0' sn4 = 'Charge Density' dim = 1 obj => this%species%spe(i) case ('jx') sn1 = 'Jx' sn2 = 'jx' sn3 = 'n_0 c' sn4 = 'J_x' dim = 1 obj => this%fields%cu3d case ('jy') sn1 = 'Jy' sn2 = 'jy' sn3 = 'n_0 c' sn4 = 'J_y' dim = 2 obj => this%fields%cu3d case ('jz') sn1 = 'Jz' sn2 = 'jz' sn3 = 'n_0 c' sn4 = 'J_z' dim = 3 obj => this%fields%cu3d end select if (this%in%found('species('//trim(s1)//').diag'//'('//trim(s2)//').slice')) then call this%in%info('species('//trim(s1)//').diag'//'('//trim(s2)//').slice',n_children=n) do ii = 1, n n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => obj allocate(this%diag(n_diag)%slice,this%diag(n_diag)%slice_pos) allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = dim if (allocated(sl)) deallocate(sl) write (s4, '(I4.4)') ii call this%in%get('species('//trim(s1)//').diag'//'('//trim(s2)//').& &slice('//trim(s4)//').(1)',sl) select case (sl) case ('yz') this%diag(n_diag)%slice = 1 call this%diag(n_diag)%file%new(& &axismin = (/aly1,alz1,alx1/),& &axismax = (/aly2,alz2,alx2/),& &axisname = (/'y ','\xi','x '/),& &axislabel = (/'y ','\xi','x '/)) case ('xz') this%diag(n_diag)%slice = 2 call this%diag(n_diag)%file%new(& &axismin = (/alx1,alz1,aly1/),& &axismax = (/alx2,alz2,aly2/),& &axisname = (/'x ','\xi','y '/),& &axislabel = (/'x ','\xi','y '/)) case ('xy') this%diag(n_diag)%slice = 3 call this%diag(n_diag)%file%new(& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &axisname = (/'x','y','z'/),& &axislabel = (/'x','y','z'/)) end select call this%in%get('species('//trim(s1)//').diag'//'('//trim(s2)//').& &slice('//trim(s4)//').(2)',this%diag(n_diag)%slice_pos) if (this%p%getidproc() == 0) then call system('mkdir -p ./Species'//trim(s1)//'/'//trim(sn1)//'_slice_'//trim(s4)//'/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &rank = 2,& &filename = './Species'//trim(s1)//'/'//trim(sn1)//'_slice_'//trim(s4)//'/',& &dataname = trim(sn2)//'_slice_'//sl,& &units = trim(sn3),& &label = trim(sn4)) end do else n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => obj allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = dim if (this%p%getidproc() == 0) then call system('mkdir -p ./Species'//trim(s1)//'/'//trim(sn1)//'/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisname = (/'x ','y ','\xi'/),& &axislabel = (/'x ','y ','\xi'/),& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &rank = 3,& &filename = './Species'//trim(s1)//'/'//trim(sn1)//'/',& &dataname = trim(sn2),& &units = trim(sn3),& &label = trim(sn4)) end if end do end if end do end do call this%in%info('field.diag',n_children=n) do i = 1, n write (s1, '(I4.4)') i call this%in%get('field.diag('//trim(s1)//').ndump',ndump) if (ndump>0) then call this%in%info('field.diag('//trim(s1)//').name',n_children=m) do j = 1, m write (s2, '(I4.4)') j if (allocated(ss)) deallocate(ss) call this%in%get('field.diag('//trim(s1)//').name('//trim(s2)//')',ss) select case (trim(ss)) case ('ex') sn1 = 'Ex' sn2 = 'ex' sn3 = 'mc\omega_p/e' sn4 = 'Electric Field' dim = 1 obj => this%fields%bexyz case ('ey') sn1 = 'Ey' sn2 = 'ey' sn3 = 'mc\omega_p/e' sn4 = 'Electric Field' dim = 2 obj => this%fields%bexyz case ('ez') sn1 = 'Ez' sn2 = 'ez' sn3 = 'mc\omega_p/e' sn4 = 'Electric Field' dim = 3 obj => this%fields%bexyz case ('bx') sn1 = 'Bx' sn2 = 'bx' sn3 = 'mc\omega_p/e' sn4 = 'Magnetic Field' dim = 1 obj => this%fields%bbxyz case ('by') sn1 = 'By' sn2 = 'by' sn3 = 'mc\omega_p/e' sn4 = 'Magnetic Field' dim = 2 obj => this%fields%bbxyz case ('bz') sn1 = 'Bz' sn2 = 'bz' sn3 = 'mc\omega_p/e' sn4 = 'Magnetic Field' dim = 3 obj => this%fields%bbxyz case ('psi') sn1 = 'Psi' sn2 = 'psi' sn3 = 'mc^2' sn4 = '\Psi' dim = 1 obj => this%fields%psi3d end select if (this%in%found('field.diag('//trim(s1)//').slice')) then call this%in%info('field.diag('//trim(s1)//').slice',n_children=l) do k = 1, l n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => obj allocate(this%diag(n_diag)%slice,this%diag(n_diag)%slice_pos) allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = dim if (allocated(sl)) deallocate(sl) write (s3, '(I4.4)') k call this%in%get('field.diag('//trim(s1)//').& &slice('//trim(s3)//').(1)',sl) select case (sl) case ('yz') this%diag(n_diag)%slice = 1 call this%diag(n_diag)%file%new(& &axismin = (/aly1,alz1,alx1/),& &axismax = (/aly2,alz2,alx2/),& &axisname = (/'y ','\xi','x '/),& &axislabel = (/'y ','\xi','x '/)) case ('xz') this%diag(n_diag)%slice = 2 call this%diag(n_diag)%file%new(& &axismin = (/alx1,alz1,aly1/),& &axismax = (/alx2,alz2,aly2/),& &axisname = (/'x ','\xi','y '/),& &axislabel = (/'x ','\xi','y '/)) case ('xy') this%diag(n_diag)%slice = 3 call this%diag(n_diag)%file%new(& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &axisname = (/'x','y','z'/),& &axislabel = (/'x','y','z'/)) end select call this%in%get('field.diag('//trim(s1)//').& &slice('//trim(s3)//').(2)',this%diag(n_diag)%slice_pos) if (this%p%getidproc() == 0) then call system('mkdir -p ./Fields/'//trim(sn1)//'_slice'//trim(s3)//'/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &rank = 2,& &filename = './Fields/'//trim(sn1)//'_slice'//trim(s3)//'/',& &dataname = trim(sn2)//'slice'//sl,& &units = trim(sn3),& &label = trim(sn4)) end do else n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => obj allocate(this%diag(n_diag)%dim) this%diag(n_diag)%dim = dim if (this%p%getidproc() == 0) then call system('mkdir -p ./Fields/'//trim(sn1)//'/') end if call this%diag(n_diag)%file%new(& &timeunits = '1 / \omega_p',& &dt = dt,& &axisname = (/'x ','y ','\xi'/),& &axislabel = (/'x ','y ','\xi'/),& &axisunits = (/'c / \omega_p','c / \omega_p','c / \omega_p'/),& &axismin = (/alx1,aly1,alz1/),& &axismax = (/alx2,aly2,alz2/),& &rank = 3,& &filename = './Fields/'//trim(sn1)//'/',& &dataname = trim(sn2),& &units = trim(sn3),& &label = trim(sn4)) end if end do end if end do if (rst) then call this%in%get('simulation.ndump_restart',ndump) do i = 1, this%nbeams write (s1,'(I4.4)') i write (s2,'(I10.10)') this%p%getidproc() n_diag = n_diag + 1 this%diag(n_diag)%df = ndump this%diag(n_diag)%obj => this%beams%beam(i) if (this%p%getidproc() == 0) then call system('mkdir -p ./RST/Beam-'//trim(s1)//'/') end if call this%diag(n_diag)%file%new(& &ty='rst',& &filename = './RST/Beam-'//trim(s1)//'/',& &dataname = 'RST-beam'//trim(s1)//'-'//trim(s2)) end do end if call MPI_BARRIER(this%p%getlworld(),ierr) call this%err%werrfl2(class//sname//' ended') end subroutine init_diag ! subroutine diag_simulation(this) implicit none class(simulation), intent(inout) :: this ! local data character(len=18), save :: sname = 'diag_simulation:' integer :: n, m, l, i, j, k, ierr, idn = 9 integer, dimension(10) :: istat real :: dt call this%err%werrfl2(class//sname//' started') call this%in%get('simulation.dt',dt) n = size(this%diag) do i = 1, n if (mod(this%tstep-1,this%diag(i)%df) == 0) then call this%diag(i)%file%new(n = this%tstep, t = this%tstep*dt) select type (obj => this%diag(i)%obj) type is (field3d) if (allocated(this%diag(i)%slice)) then this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wr(this%diag(i)%file,this%diag(i)%dim,this%diag(i)%slice,this%diag(i)%slice_pos,& &this%tag(1),this%tag(1),this%id(idn+i)) else this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wr(this%diag(i)%file,this%diag(i)%dim,this%tag(1),this%tag(1),this%id(idn+i)) end if type is (beam3d) if (allocated(this%diag(i)%psample)) then this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wr(this%diag(i)%file,this%diag(i)%psample,(/this%dex,this%dex,this%dxi/),& &this%tag(1),this%tag(1),this%id(idn+i)) else if (allocated(this%diag(i)%slice)) then this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wrq(this%diag(i)%file,this%diag(i)%slice,this%diag(i)%slice_pos,& &this%tag(1),this%tag(1),this%id(idn+i)) else if (allocated(this%diag(i)%dim)) then this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wrq(this%diag(i)%file,this%tag(1),this%tag(1),this%id(idn+i)) else call obj%wrst(this%diag(i)%file) end if type is (species2d) if (allocated(this%diag(i)%slice)) then this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wrq(this%diag(i)%file,this%diag(i)%slice,this%diag(i)%slice_pos,& &this%tag(1),this%tag(1),this%id(idn+i)) else this%tag(1) = ntag() call MPI_WAIT(this%id(idn+i),istat,ierr) call obj%wrq(this%diag(i)%file,this%tag(1),this%tag(1),this%id(idn+i)) end if end select end if end do call this%err%werrfl2(class//sname//' ended') end subroutine diag_simulation ! function ntag() implicit none integer, save :: tag = 0 integer :: ntag ntag = tag tag = tag + 1 end function ntag ! end module simulation_class