diff --git a/input_file/analytic_function/qpinput.json b/input_file/analytic_function/qpinput.json new file mode 100644 index 0000000..d367352 --- /dev/null +++ b/input_file/analytic_function/qpinput.json @@ -0,0 +1,189 @@ +{ + "simulation": { + "algorithm": "standard", + "nodes": [ + 8, + 2 + ], + "indx": 9, + "indy": 9, + "indz": 9, + "box": { + "x": [ + -30, + 30 + ], + "y": [ + -30, + 30 + ], + "z": [ + 0, + 9.45 + ] + }, + "boundary": "conducting", + "interpolation": "linear", + "n0": 1e+18, + "time": 101, + "dt": 10.0, + "nbeams": 1, + "nspecies": 1, + "dump_restart": false, + "ndump_restart": 1, + "read_restart": false, + "restart_timestep": 0, + "iter": 1, + "verbose": 0 + }, + "beam": [ + { + "evolution": true, + "profile": 13, + "ppc": [ + 2, + 2, + 2 + ], + "density" : 1.0, + "npmax": 234217728, + "q": -1.0, + "m": 1.0, + "gamma": 97847, + "quiet_start": true, + "center": [ + 0.0, + 0.0, + 0.0 + ], + "math_func" : "exp(-(z-3)^2/1.0) * exp(-(x^2 + y^2)/1.0)", + "sigma_v": [ + 0.0, + 0.0, + 0.0 + ], + "range_x": [ + -3.0, + 3.0 + ], + "range_y": [ + -3.0, + 3.0 + ], + "range_z" : [ + 1, + 5 + ], + "centroid_x": [ + 0.0, + 0.0, + 0.0 + ], + "centroid_y": [ + 0.0, + 0.0, + 0.0 + ], + "diag": [ + { + "name": [ + "charge" + ], + "ndump": 0 + }, + { + "name": [ + "charge" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 257 + ], + [ + "yz", + 257 + ] + ] + }, + { + "name": [ + "raw" + ], + "ndump": 1, + "sample": 100 + } + ] + } + ], + "species": [ + { + "profile": 13, + "ppc": [ + 2, + 2 + ], + "q": -1.0, + "m": 1.0, + !"push_type" : "standard", + "push_type" : "robust", + "math_func" : "s^2/(s + 10.0)^2 * exp(-(x^2 + y^2)/50.0)", + "density": 1.0, + "npmax": 3200000, + "diag": [ + { + "name": [ + "charge" + ], + "ndump": 0 + }, + { + "name": [ + "charge","jx","jy","jz" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 256 + ], + [ + "yz", + 256 + ] + ] + } + ] + } + ], + "field": { + "diag": [ + { + "name": [ + "ex", + "ez" + ], + "ndump": 0 + }, + { + "name": [ + "ez", "ex", "ey", "bx", "by", "bz","psi" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 256 + ], + [ + "yz", + 256 + ] + ] + } + ] + } +} + + diff --git a/input_file/robust/README.md b/input_file/robust/README.md new file mode 100644 index 0000000..5c2ed63 --- /dev/null +++ b/input_file/robust/README.md @@ -0,0 +1,2 @@ +# QuickPIC Input Decks +A sample input deck using the robust pusher. diff --git a/input_file/robust/qpinput.json b/input_file/robust/qpinput.json new file mode 100644 index 0000000..e9fa09d --- /dev/null +++ b/input_file/robust/qpinput.json @@ -0,0 +1,180 @@ +{ + "simulation": { + "algorithm": "standard", + "nodes": [ + 1, + 1 + ], + "indx": 9, + "indy": 9, + "indz": 9, + "box": { + "x": [ + -30, + 30 + ], + "y": [ + -30, + 30 + ], + "z": [ + 0, + 9.45 + ] + }, + "boundary": "conducting", + "interpolation": "linear", + "n0": 1e+18, + "time": 11, + "dt": 10.0, + "nbeams": 1, + "nspecies": 1, + "dump_restart": false, + "ndump_restart": 1, + "read_restart": false, + "restart_timestep": 0, + "iter": 1, + "verbose": 0 + }, + "beam": [ + { + "evolution": true, + "profile": 0, + "np": [ + 128, + 128, + 128 + ], + "npmax": 32000000, + "q": -1.0, + "m": 1.0, + "gamma": 97847, + "peak_density": 400.0, + "quiet_start": true, + "center": [ + 0.0, + 0.0, + 1.0 + ], + "sigma": [ + 0.18, + 0.18, + 0.18818 + ], + "sigma_v": [ + 0.0, + 0.0, + 0.0 + ], + "centroid_x": [ + 0.0, + 0.0, + 0.0 + ], + "centroid_y": [ + 0.0, + 0.0, + 0.0 + ], + "diag": [ + { + "name": [ + "charge" + ], + "ndump": 0 + }, + { + "name": [ + "charge" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 256 + ], + [ + "yz", + 256 + ] + ] + }, + { + "name": [ + "raw" + ], + "ndump": 1, + "sample": 100 + } + ] + } + ], + "species": [ + { + "profile": 0, + "ppc": [ + 2, + 2 + ], + "q": -1.0, + "m": 1.0, + !"push_type" : "standard", + "push_type" : "robust", + "longitudinal_profile" : "uniform", + "density": 1.0, + "npmax": 32000000, + "diag": [ + { + "name": [ + "charge" + ], + "ndump": 0 + }, + { + "name": [ + "charge","jx","jy","jz" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 256 + ], + [ + "yz", + 256 + ] + ] + } + ] + } + ], + "field": { + "diag": [ + { + "name": [ + "ex", + "ez" + ], + "ndump": 0 + }, + { + "name": [ + "ez", "ex", "ey", "bx", "by", "bz","psi" + ], + "ndump": 1, + "slice": [ + [ + "xz", + 256 + ], + [ + "yz", + 256 + ] + ] + } + ] + } +} + diff --git a/source/fdist2d_class.f03 b/source/fdist2d_class.f03 index 4fec744..441e224 100644 --- a/source/fdist2d_class.f03 +++ b/source/fdist2d_class.f03 @@ -8,12 +8,12 @@ module fdist2d_class use spect2d_class use ufield2d_class use input_class - + use m_fparser implicit none private - public :: fdist2d, fdist2d_000, fdist2d_010, fdist2d_011, fdist2d_012 + public :: fdist2d, fdist2d_000, fdist2d_010, fdist2d_011, fdist2d_012, fdist2d_013 type, abstract :: fdist2d @@ -126,6 +126,22 @@ end subroutine ab_init_fdist2d procedure, private :: dist2d => dist2d_012 end type fdist2d_012 + + type, extends(fdist2d) :: fdist2d_013 +! Transeversely uniform profile with uniform or piecewise longitudinal profile + private +! xppc, yppc = particle per cell in x and y directions + integer :: xppc, yppc + real :: qm, den + real :: cx, cy, dx, dy + ! analytic density math function + type(t_fparser), pointer :: math_func => null() + + contains + procedure, private :: init_fdist2d => init_fdist2d_013 + procedure, private :: dist2d => dist2d_013 + + end type fdist2d_013 ! character(len=10), save :: class = 'fdist2d:' character(len=128), save :: erstr @@ -949,4 +965,176 @@ subroutine dist2d_012(this,part2d,npp,fd,s) end subroutine dist2d_012 ! + + subroutine init_fdist2d_013(this,input,i) + + implicit none + + class(fdist2d_013), intent(inout) :: this + type(input_json), intent(inout), pointer :: input + integer, intent(in) :: i +! local data + integer :: npf,xppc,yppc,npmax,indx,indy + real :: qm, den + real :: cx, cy + real :: min, max + real :: alx, aly, dx, dy + character(len=20) :: sn,s1 + character(len=18), save :: sname = 'init_fdist2d_013:' + character(len=:), allocatable :: read_str + integer :: ierr + this%sp => input%sp + this%err => input%err + this%p => input%pp + + call this%err%werrfl2(class//sname//' started') + write (sn,'(I3.3)') i + s1 = 'species('//trim(sn)//')' + call input%get('simulation.indx',indx) + call input%get('simulation.indy',indy) + call input%get('simulation.box.x(1)',min) + call input%get('simulation.box.x(2)',max) + cx = - min + alx = (max-min) + dx=alx/real(2**indx) + call input%get('simulation.box.y(1)',min) + call input%get('simulation.box.y(2)',max) + cy = - min + aly = (max-min) + dy=aly/real(2**indy) + call input%get(trim(s1)//'.profile',npf) + call input%get(trim(s1)//'.ppc(1)',xppc) + call input%get(trim(s1)//'.ppc(2)',yppc) + call input%get(trim(s1)//'.q',qm) + call input%get(trim(s1)//'.density',den) + npmax = xppc*yppc*(2**indx)*(2**indy)/this%p%getlnvp()*4 + + this%dx = dx + this%dy = dy + this%npf = npf + this%xppc = xppc + this%yppc = yppc + this%qm = qm + this%den = den + this%npmax = npmax + this%cx = cx/dx + this%cy = cy/dy + + if ( .not. associated( this%math_func ) ) then + allocate( t_fparser :: this%math_func ) + endif + call input%get( trim(s1) // '.math_func', read_str ) + call setup(this%math_func, trim(read_str), (/'x','y','s'/), ierr) + + call this%err%werrfl2(class//sname//' ended') + + + end subroutine init_fdist2d_013 +! + subroutine dist2d_013(this,part2d,npp,fd,s) + implicit none + class(fdist2d_013), intent(inout) :: this + real, dimension(:,:), pointer, intent(inout) :: part2d + integer, intent(inout) :: npp + class(ufield2d), intent(in), pointer :: fd + real, intent(in) :: s +! local data + character(len=18), save :: sname = 'dist2d_013:' + real, dimension(:,:), pointer :: pt => null() + integer :: nps, nx, ny, noff, xppc, yppc, i, j + integer :: ix, iy + real :: qm, x, y + real :: cx, cy + real :: den_temp, rr + integer :: prof_l, ii + real(p_k_fparse), dimension(3) :: fparser_arr + + call this%err%werrfl2(class//sname//' started') + + nx = fd%getnd1p(); ny = fd%getnd2p(); noff = fd%getnoff() + xppc = this%xppc; yppc = this%yppc + den_temp = 1.0 + qm = den_temp*this%den*this%qm/abs(this%qm)/real(xppc)/real(yppc) + cx = this%cx; cy = this%cy + nps = 1 + pt => part2d + fparser_arr(3) = s +! initialize the particle positions + if (noff < ny) then + do i=2, nx-1 + do j=2, ny + do ix = 0, xppc-1 + do iy=0, yppc-1 + x = (ix + 0.5)/xppc + i - 1 + y = (iy + 0.5)/yppc + j - 1 + noff + fparser_arr(1) = (x-cx) * this%dx + fparser_arr(2) = (y-cy) * this%dy + den_temp = eval(this%math_func, fparser_arr) + pt(1,nps) = x + pt(2,nps) = y + pt(3,nps) = 0.0 + pt(4,nps) = 0.0 + pt(5,nps) = 0.0 + pt(6,nps) = 1.0 + pt(7,nps) = 1.0 + pt(8,nps) = qm*den_temp + nps = nps + 1 + enddo + enddo + enddo + enddo + else if (noff > (nx-ny-1)) then + do i=2, nx-1 + do j=1, ny-1 + do ix = 0, xppc-1 + do iy=0, yppc-1 + x = (ix + 0.5)/xppc + i - 1 + y = (iy + 0.5)/yppc + j - 1 + noff + fparser_arr(1) = (x-cx) * this%dx + fparser_arr(2) = (y-cy) * this%dy + den_temp = eval(this%math_func, fparser_arr) + pt(1,nps) = x + pt(2,nps) = y + pt(3,nps) = 0.0 + pt(4,nps) = 0.0 + pt(5,nps) = 0.0 + pt(6,nps) = 1.0 + pt(7,nps) = 1.0 + pt(8,nps) = qm*den_temp + nps = nps + 1 + enddo + enddo + enddo + enddo + else + do i=2, nx-1 + do j=1, ny + do ix = 0, xppc-1 + do iy=0, yppc-1 + x = (ix + 0.5)/xppc + i - 1 + y = (iy + 0.5)/yppc + j - 1 + noff + fparser_arr(1) = (x-cx) * this%dx + fparser_arr(2) = (y-cy) * this%dy + den_temp = eval(this%math_func, fparser_arr) + pt(1,nps) = x + pt(2,nps) = y + pt(3,nps) = 0.0 + pt(4,nps) = 0.0 + pt(5,nps) = 0.0 + pt(6,nps) = 1.0 + pt(7,nps) = 1.0 + pt(8,nps) = qm*den_temp + nps = nps + 1 + enddo + enddo + enddo + enddo + endif + + npp = nps - 1 + + call this%err%werrfl2(class//sname//' ended') + + end subroutine dist2d_013 + end module fdist2d_class \ No newline at end of file diff --git a/source/fdist3d_class.f03 b/source/fdist3d_class.f03 index 54c5327..c984f51 100644 --- a/source/fdist3d_class.f03 +++ b/source/fdist3d_class.f03 @@ -9,6 +9,7 @@ module fdist3d_class use ufield3d_class use part3d_lib use input_class + use m_fparser implicit none @@ -16,6 +17,7 @@ module fdist3d_class public :: fdist3d, fdist3d_000, fdist3d_001, fdist3d_002, fdist3d_100 public :: fdist3d_003 + public :: fdist3d_013 type, abstract :: fdist3d @@ -131,6 +133,22 @@ end subroutine ab_init_fdist3d end type fdist3d_003 ! + type, extends(fdist3d) :: fdist3d_013 +! math func profile + private + + integer :: xppc, yppc, zppc + real :: qm, bcx, bcy, bcz, dx, dy, dz + real, dimension(3) :: lb, ub + real :: gamma, np + type(t_fparser), pointer :: math_func => null() + + contains + procedure, private :: init_fdist3d => init_fdist3d_013 + procedure, private :: dist3d => dist3d_013 + + end type fdist3d_013 + type, extends(fdist3d) :: fdist3d_100 ! Ring profile private @@ -972,6 +990,202 @@ subroutine dist3d_003(this,part3d,npp,fd) call this%err%werrfl2(class//sname//' ended') end subroutine dist3d_003 + + subroutine init_fdist3d_013(this,input,i) + + implicit none + + class(fdist3d_013), intent(inout) :: this + type(input_json), intent(inout), pointer :: input + integer, intent(in) :: i +! local data + integer :: xppc,yppc,zppc,npmax,npf + real :: qm, bcx, bcy, bcz + real, dimension(3) :: lb, ub + real :: gamma, np + real :: min, max, n0 + real :: alx, aly, alz, dx, dy, dz + integer :: indx, indy, indz, ierr + character(len=20) :: sn,s1 + + character(len=18), save :: sname = 'init_fdist3d_013:' + character(len=:), allocatable :: read_str + + this%sp => input%sp + this%err => input%err + this%p => input%pp + + call this%err%werrfl2(class//sname//' started') + + write (sn,'(I3.3)') i + s1 = 'beam('//trim(sn)//')' + + call input%get('simulation.n0',n0) + call input%get('simulation.indx',indx) + call input%get('simulation.indy',indy) + call input%get('simulation.indz',indz) + + call input%get('simulation.box.x(1)',min) + call input%get('simulation.box.x(2)',max) + call input%get(trim(s1)//'.center(1)',bcx) + call input%get(trim(s1)//'.range_x(1)',lb(1)) + call input%get(trim(s1)//'.range_x(2)',ub(1)) + lb(1) = lb(1) - min + ub(1) = ub(1) - min + bcx = bcx - min + alx = (max-min) + dx=alx/real(2**indx) + call input%get('simulation.box.y(1)',min) + call input%get('simulation.box.y(2)',max) + call input%get(trim(s1)//'.center(2)',bcy) + call input%get(trim(s1)//'.range_y(1)',lb(2)) + call input%get(trim(s1)//'.range_y(2)',ub(2)) + lb(2) = lb(2) - min + ub(2) = ub(2) - min + bcy = bcy -min + aly = (max-min) + dy=aly/real(2**indy) + call input%get('simulation.box.z(1)',min) + call input%get('simulation.box.z(2)',max) + call input%get(trim(s1)//'.center(3)',bcz) + call input%get(trim(s1)//'.range_z(1)',lb(3)) + call input%get(trim(s1)//'.range_z(2)',ub(3)) + lb(3) = lb(3) - min + ub(3) = ub(3) - min + bcz = bcz -min + alz = (max-min) + dz=alz/real(2**indz) + + call input%get(trim(s1)//'.profile',npf) + call input%get(trim(s1)//'.ppc(1)',xppc) + call input%get(trim(s1)//'.ppc(2)',yppc) + call input%get(trim(s1)//'.ppc(3)',zppc) + call input%get(trim(s1)//'.q',qm) + call input%get(trim(s1)//'.density',np) + call input%get(trim(s1)//'.gamma',gamma) + call input%get(trim(s1)//'.npmax',npmax) + + this%npf = npf + this%npmax = npmax + this%xppc = xppc + this%yppc = yppc + this%zppc = zppc + this%bcx = bcx/dx + this%bcy = bcy/dy + this%bcz = bcz/dz + this%lb(1) = lb(1)/dx + this%lb(2) = lb(2)/dy + this%lb(3) = lb(3)/dz + this%ub(1) = ub(1)/dx + this%ub(2) = ub(2)/dy + this%ub(3) = ub(3)/dz + this%dx = dx + this%dy = dy + this%dz = dz + this%gamma = gamma + this%np = np + + this%qm = qm + + if ( .not. associated( this%math_func ) ) then + allocate( t_fparser :: this%math_func ) + endif + call input%get( trim(s1) // '.math_func', read_str ) + call setup(this%math_func, trim(read_str), (/'x','y','z'/), ierr) + + call this%err%werrfl2(class//sname//' ended') + + end subroutine init_fdist3d_013 +! + subroutine dist3d_013(this,part3d,npp,fd) + + implicit none + + class(fdist3d_013), intent(inout) :: this + real, dimension(:,:), pointer, intent(inout) :: part3d + integer, intent(inout) :: npp + class(ufield3d), intent(in), pointer :: fd +! local data1 + real, dimension(:,:), pointer :: pt => null() + integer :: xppc,yppc,zppc + real :: qm, bcx, bcy, bcz + real, dimension(3) :: lb, ub + real :: r1, r2, tx, ty, tz, a1, a2, a3, a4 + real :: gamma, np, den_temp + integer, dimension(2) :: noff + integer :: nps + integer :: i, j, k, ix, iy, iz + integer :: idimp, npmax, ierr + character(len=18), save :: sname = 'dist3d_013:' + real(p_k_fparse), dimension(3) :: fparser_arr + + call this%err%werrfl2(class//sname//' started') + + nps = 1 + ierr = 0 + xppc = this%xppc; yppc = this%yppc; zppc = this%zppc; + noff = fd%getnoff() + lb(1) = max(this%lb(1),0.0) + ub(1) = min(this%ub(1),real(fd%getnd1()-1)) + + if (this%lb(2) > real(noff(1)+fd%getnd2p()) .or. & + &this%ub(2) < real(noff(1))) then + npp = 0 + return + else + lb(2) = max(this%lb(2),real(noff(1))) + ub(2) = min(this%ub(2),real(noff(1)+fd%getnd2p()-1)) + end if + if (this%lb(3) > real(noff(2)+fd%getnd3p()) .or. & + &this%ub(3) < real(noff(2))) then + npp = 0 + return + else + lb(3) = max(this%lb(3),real(noff(2))) + ub(3) = min(this%ub(3),real(noff(2)+fd%getnd3p()-1)) + end if + pt => part3d + np = this%np + gamma = this%gamma + bcx = this%bcx; bcy = this%bcy; bcz = this%bcz + qm = this%qm + den_temp = 1.0 + qm = this%np*this%qm/abs(this%qm)/real(xppc)/real(yppc)/real(zppc) + + do i=int(lb(1)), int(ub(1)) + do j=int(lb(2)), int(ub(2)) + do k = int(lb(3)), int(ub(3)) + do ix = 0, xppc-1 + do iy= 0, yppc-1 + do iz= 0, zppc-1 + tx = (ix + 0.5)/xppc + i + ty = (iy + 0.5)/yppc + j + tz = (iz + 0.5)/zppc + k + fparser_arr(1) = (tx-bcx) * this%dx + fparser_arr(2) = (ty-bcy) * this%dy + fparser_arr(3) = (tz-bcz) * this%dz + den_temp = eval(this%math_func, fparser_arr) + pt(1,nps) = tx + pt(2,nps) = ty + pt(3,nps) = tz + pt(4,nps) = 0.0 + pt(5,nps) = 0.0 + pt(6,nps) = gamma + pt(7,nps) = qm*den_temp + nps = nps + 1 + enddo + enddo + enddo + enddo + enddo + enddo + + npp = nps - 1 + + call this%err%werrfl2(class//sname//' ended') + + end subroutine dist3d_013 + ! subroutine init_fdist3d_100(this,input,i) diff --git a/source/fparser.f03 b/source/fparser.f03 new file mode 100644 index 0000000..1709010 --- /dev/null +++ b/source/fparser.f03 @@ -0,0 +1,1609 @@ +! +! Function parsing class +! + +! There is a known issue with variable names: if a variable name ends in digit+'e' (or 'd') it +! i.e. 'a2e' may be interpreted as a numeric constant and evaluating something like var +/- number +! i.e. 'a2e - 1' may give the wrong result. + +module m_fparser + + ! use sysutil_module + + implicit none + +! restrict access to things explicitly declared public + private + + ! precision to use for the calculations + !integer, parameter, public :: p_k_fparse = p_single + integer, parameter, public :: p_k_fparse = 8 + + character, parameter :: p_tab = achar(9) + character, parameter ::p_space = ' ' + + ! predefined sizes + integer, parameter :: p_max_expr_len = 1024 ! maximum length an expression can have + integer, parameter :: MAX_FUNC_LEN = 16 ! maximum length a single function name can have + integer, parameter :: MAX_VAR_LEN = 128 ! maximum length a variable name can have + integer, parameter :: MAX_VARS = 16 ! maximum number of variables + integer, parameter :: MAX_NESTED_FUNC = 16 ! maximum number of nested functions allowed + integer, parameter :: MAX_OP_LEN = 2 ! maximum length an operator can have + + ! predefined bytecode and stack sizes + integer, parameter :: MAX_CODE_SIZE = 256 ! code will be limited to 256 expressions + integer, parameter :: MAX_STACK_SIZE = 256 ! maximum number of values in the calculation stack + integer, parameter :: MAX_DATA_SIZE = 64 ! maximum number of values in the data stack + + ! error codes + integer, parameter :: err_eof = -1 + integer, parameter :: err_mult_op = -2 + integer, parameter :: err_inv_num = -3 + integer, parameter :: err_parent = -4 + integer, parameter :: err_syntax = -5 + integer, parameter :: err_fun_par = -6 + integer, parameter :: err_param = -7 + + + ! List of valid operators: + ! - operators are specified in groups of operations with the same priority, separated by + ! blank (' '), elements. + ! - operator groups are ordered with increasing priority. + ! - Note that if two operators within the same group begin with the same character and + ! have different lengths (e.g. '<=' and '<'), the longest one MUST be specified first + ! because of the implementation of the scan_ops function + ! + ! NOTE: Operators must not have alphanumeric characters or dots + + character(len = *), dimension(19), parameter :: oper = & + (/ '||', ' ', & + '&&', ' ', & + '==', '!=', ' ', & + '<=', '>=', '< ', '> ', ' ', & + '+ ', '- ', ' ', & + '* ', '/ ', ' ', & + '^ ' /) + + ! NOTE2: operator ^ is used for a ^ int(b) i.e. b is converted to an integer before + ! evaluating. If you want b to not be an integer use the function pow(a,b) instead but + ! keep in mind that if a < 0 pow(a,b) returns NaN + + ! Operator op codes + ! These must be in sequence with the same order as the operator list array + + integer, parameter :: op_or = 1, & + op_and = 2, & + op_eq = 3, op_ne = 4, & + op_le = 5, op_ge = 6, op_lt = 7, op_gt = 8, & + op_add = 9, op_sub = 10, & + op_mult = 11, op_div = 12, & + op_pow = 13 + + ! Total number of operators + integer, parameter :: n_oper = 13 + + integer, parameter :: bcode_data = 0 + + ! valid functions + integer, parameter :: n_functions = 30 + + type :: t_math_func + character (len = MAX_FUNC_LEN) :: name + integer :: n_params + end type t_math_func + + ! note that if two functions begin with the same string + ! the longest one MUST be specified first because of the + ! implementation of the function if_function, e.g. log10 and log, + ! or atan2 and atan + + type (t_math_func), dimension(n_functions) , parameter :: functions = & + (/ t_math_func('abs ', 1) , & ! 1 + t_math_func('sinh ', 1) , & ! 2 + t_math_func('sin ', 1) , & ! 3 + t_math_func('cosh ', 1) , & ! 4 + t_math_func('cos ', 1) , & ! 5 + t_math_func('tanh ', 1) , & ! 6 + t_math_func('tan ', 1) , & ! 7 + t_math_func('exp ', 1) , & ! 8 + t_math_func('log10 ', 1) , & ! 9 + t_math_func('log ', 1) , & ! 10 + t_math_func('asin ', 1) , & ! 11 + t_math_func('acos ', 1) , & ! 12 + t_math_func('atan2 ', 2) , & ! 13 + t_math_func('atan ', 1) , & ! 14 + t_math_func('sqrt ', 1) , & ! 15 + t_math_func('not ', 1) , & ! 16 + t_math_func('neg ', 1) , & ! 17 + t_math_func('if ', 3) , & ! 18 + t_math_func('pow ', 2) , & ! 19 + t_math_func('int ', 1) , & ! 20 + t_math_func('nint ', 1) , & ! 21 + t_math_func('ceiling ', 1) , & ! 22 + t_math_func('floor ', 1) , & ! 23 + t_math_func('modulo ', 2) , & ! 24 + t_math_func('rect ', 1) , & ! 25 + t_math_func('step ', 1) , & ! 26 + t_math_func('min3 ', 3) , & ! 27 + t_math_func('min ', 2) , & ! 28 + t_math_func('max3 ', 3) , & ! 29 + t_math_func('max ', 2) /) ! 30 + + integer, parameter :: func_abs = 1, & + func_sinh = 2, & + func_sin = 3, & + func_cosh = 4, & + func_cos = 5, & + func_tanh = 6, & + func_tan = 7, & + func_exp = 8, & + func_log10 = 9, & + func_log = 10, & + func_asin = 11, & + func_acos = 12, & + func_atan2 = 13, & + func_atan = 14, & + func_sqrt = 15, & + func_not = 16, & + func_neg = 17, & + func_if = 18, & + func_pow = 19, & + func_int = 20, & + func_nint = 21, & + func_ceiling = 22, & + func_floor = 23, & + func_modulo = 24, & + func_rect = 25, & + func_step = 26, & + func_min3 = 27, & + func_min = 28, & + func_max3 = 29, & + func_max = 30 + + + ! beginning of variable space code + integer, parameter :: bcode_var = n_oper + n_functions + + ! class definition + type :: t_fparser + + integer :: code_size ! actual size of compiled code + integer :: stack_size ! actual size of calculation stack + integer :: data_size ! actual size of data stack + + integer :: stack_ptr ! pointer to the current value on the stack + ! used at compile time to determine required stack size + ! since we are not allocating the stack dynamically it + ! serves no purpose + + ! compiled code + integer, dimension(MAX_CODE_SIZE) :: bytecode + + ! stack for temporary values during calculation + real(p_k_fparse), dimension(MAX_STACK_SIZE):: stack + + ! internal stack for compiled number constants + real(p_k_fparse), dimension(MAX_DATA_SIZE):: data_stack + + ! Expression text + character (len = p_max_expr_len) :: func + + ! Number of variables + integer :: nvars + + ! Variable names + character (len = MAX_VAR_LEN), dimension(MAX_VARS) :: var_list + + end type t_fparser + + interface setup + module procedure setup_fparser + end interface + + interface cleanup + module procedure cleanup_fparser + end interface + + interface eval + module procedure eval_fparser + end interface + + interface functext + module procedure functext_parser + end interface + +! declare things that should be public + + public :: t_fparser, setup, cleanup, eval, p_max_expr_len, functext + + contains + + function functext_parser( this ) + implicit none + type( t_fparser ), intent(in) :: this + character( len = len_trim( this%func ) ) :: functext_parser + + functext_parser = trim( this%func ) + end function functext_parser + + + + +!----------------------------------------------------------------------------------------- +! class constructor +!----------------------------------------------------------------------------------------- + subroutine setup_fparser( this, str_expr, var_list, ierr ) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + character ( len = * ), intent (in) :: str_expr + character ( len = * ), dimension(:), intent (in) :: var_list + integer, intent (out) :: ierr + + ! executable statements + + + ierr = 0 + + ! initialize local data + + this%func = remove_spaces(str_expr) + this%nvars = size(var_list) + this%var_list(1:this%nvars) = var_list + + ! not necessary + if (this%nvars < MAX_VARS) this%var_list(this%nvars+1:) = '~' + + ! check expression syntax + + call check_syntax_fparser( this , ierr ) + if (ierr < 0) return + + ! compile the function + + call compile_fparser( this ) + + end subroutine setup_fparser +!----------------------------------------------------------------------------------------- + +subroutine cleanup_fparser( this ) + + implicit none + + type( t_fparser ), intent(inout) :: this + + this%func = '' + this%nvars = 0 + this%var_list = '' + + this%code_size = 0 + +end subroutine cleanup_fparser + +!----------------------------------------------------------------------------------------- +! process bytecode and evaluate function +!----------------------------------------------------------------------------------------- + function eval_fparser( this, var_values ) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + real( p_k_fparse ) :: eval_fparser + type( t_fparser ), intent(inout) :: this + real( p_k_fparse ), dimension(:), intent(in) :: var_values + + ! local variables + + integer :: ip, & ! instruction pointer + sp, & ! stack pointer + dp ! data pointer + + integer :: exponent ! workarround for bug in a**b + + ! executable statements + + sp = 0; dp = 1 + + do ip = 1, this%code_size + select case (this%bytecode(ip)) + + case (bcode_data) ! compiled value, grab it from the data stack + sp = sp + 1; this%stack(sp) = this%data_stack(dp); dp = dp + 1 + + case ( op_eq ) + if ( this%stack(sp-1) == this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_ne ) + if ( this%stack(sp-1) /= this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_ge ) + if ( this%stack(sp-1) >= this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_le ) + if ( this%stack(sp-1) <= this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_gt ) + if ( this%stack(sp-1) > this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_lt ) + if ( this%stack(sp-1) < this%stack(sp) ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_and ) + if ( this%stack(sp-1) /= 0.0 .and. this%stack(sp) /= 0.0 ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_or ) + if ( this%stack(sp-1) /= 0.0 .or. this%stack(sp) /= 0.0 ) then + this%stack(sp-1) = 1.0 + else + this%stack(sp-1) = 0.0 + endif + sp = sp-1 + + case ( op_add ) + this%stack(sp-1) = this%stack(sp-1) + this%stack(sp); sp = sp -1 + + case ( op_sub ) + this%stack(sp-1) = this%stack(sp-1) - this%stack(sp); sp = sp -1 + + case ( op_mult ) + this%stack(sp-1) = this%stack(sp-1) * this%stack(sp); sp = sp -1 + + case ( op_div ) + this%stack(sp-1) = this%stack(sp-1) / this%stack(sp); sp = sp -1 + + case ( op_pow ) + exponent = int(this%stack(sp)) + sp = sp -1 + this%stack(sp) = this%stack(sp) ** exponent + + case ( n_oper + func_abs ) + this%stack(sp) = abs( this%stack(sp) ) + + case ( n_oper + func_sin ) + this%stack(sp) = sin( this%stack(sp) ) + + case ( n_oper + func_sinh ) + this%stack(sp) = sinh( this%stack(sp) ) + + case ( n_oper + func_cos ) + this%stack(sp) = cos( this%stack(sp) ) + + case ( n_oper + func_cosh ) + this%stack(sp) = cosh( this%stack(sp) ) + + case ( n_oper + func_tan ) + this%stack(sp) = tan( this%stack(sp) ) + + case ( n_oper + func_tanh ) + this%stack(sp) = tanh( this%stack(sp) ) + + case ( n_oper + func_exp ) + this%stack(sp) = exp( this%stack(sp) ) + + case ( n_oper + func_log10 ) + this%stack(sp) = log10( this%stack(sp) ) + + case ( n_oper + func_log ) + this%stack(sp) = log( this%stack(sp) ) + + case ( n_oper + func_asin ) + this%stack(sp) = asin( this%stack(sp) ) + + case ( n_oper + func_acos ) + this%stack(sp) = acos( this%stack(sp) ) + + case ( n_oper + func_atan2 ) + this%stack(sp-1) = atan2(this%stack(sp-1),this%stack(sp)); sp = sp -1 + + case ( n_oper + func_atan ) + this%stack(sp) = atan( this%stack(sp) ) + + case ( n_oper + func_sqrt ) + this%stack(sp) = sqrt( this%stack(sp) ) + + case ( n_oper + func_not ) + if (this%stack(sp) == 0.0) then + this%stack(sp) = 1.0 + else + this%stack(sp) = 0.0 + endif + case ( n_oper + func_neg ) + this%stack(sp) = - this%stack(sp) + + case ( n_oper + func_if ) + if (this%stack(sp-2) /= 0.0) then + this%stack(sp-2) = this%stack(sp-1) + else + this%stack(sp-2) = this%stack(sp) + endif + sp = sp - 2 + + case ( n_oper + func_pow ) + this%stack(sp-1) = this%stack(sp-1) ** this%stack(sp); sp = sp -1 + + case ( n_oper + func_int ) + this%stack(sp) = int( this%stack(sp) ) + + case ( n_oper + func_nint ) + this%stack(sp) = nint( this%stack(sp) ) + + case ( n_oper + func_ceiling ) + this%stack(sp) = ceiling( this%stack(sp) ) + + case ( n_oper + func_floor ) + this%stack(sp) = floor( this%stack(sp) ) + + case ( n_oper + func_modulo ) + this%stack(sp-1) = modulo(this%stack(sp-1),this%stack(sp)); sp = sp -1 + + case ( n_oper + func_rect ) + if ( abs(this%stack(sp)) <= 0.5 ) then + this%stack(sp) = 1.0 + else + this%stack(sp) = 0.0 + endif + + case ( n_oper + func_step ) + if ( this%stack(sp) >= 0.0 ) then + this%stack(sp) = 1.0 + else + this%stack(sp) = 0.0 + endif + + case ( n_oper + func_min3 ) + this%stack(sp-2) = & + min(this%stack(sp-2),this%stack(sp-1),this%stack(sp)) + sp = sp -2 + + case ( n_oper + func_min ) + this%stack(sp-1) = & + min(this%stack(sp-1),this%stack(sp)); sp = sp -1 + + case ( n_oper + func_max3 ) + this%stack(sp-2) = & + max(this%stack(sp-2),this%stack(sp-1),this%stack(sp)) + sp = sp -2 + + case ( n_oper + func_max ) + this%stack(sp-1) = & + max(this%stack(sp-1),this%stack(sp)); sp = sp -1 + + case default ! this is only reached for variables + sp = sp + 1 + this%stack(sp) = var_values(this%bytecode(ip) - bcode_var ) + end select + enddo + + + eval_fparser = this%stack(1) + + end function eval_fparser +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! compile the supplied function +!----------------------------------------------------------------------------------------- + subroutine compile_fparser( this ) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + + ! local variables + + this%code_size = 0 + this%stack_size = 0 + this%data_size = 0 + + this%stack_ptr = 0 + + call compile_substr_fparser( this, trim(this%func) ) + + ! (*debug*) + !print *, "compiled code size = ", this%code_size + !print *, "required stack space = ", this%stack_size + !print *, "required data size = ", this%data_size + !print *, "code = ", this%bytecode(1:this%code_size) + + end subroutine compile_fparser +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! recursively compile substring +!----------------------------------------------------------------------------------------- +recursive subroutine compile_substr_fparser( this, sub_str ) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + character ( len = * ), intent (in) :: sub_str + + ! local variables + + integer :: sub_str_len ! length of the sub_str begin processed + integer :: fidx, f_len, n_params ! function index, function name length, number of params + integer :: varidx, var_len + + integer :: i, j, & ! string position index + p, & ! parameter count + k ! parenthesis count + + integer :: op_code, & ! offset for calculating op code + n_op, & ! number of operators in group + op_len, & ! length of operator found + list_idx, & ! position in the operator list + op_idx ! position of the operator in the operator group + + integer :: ierr + + ! print *, '>', trim(sub_str), '<' + + sub_str_len = len(sub_str) + + ! check for special cases of substring + + if (sub_str(1:1) == '+') then ! case 1: str = '+...' + call compile_substr_fparser( this, sub_str(2:) ) ! no compilation required + return + + elseif (if_enclosed_parent(sub_str)) then ! case 2: str = '(...)' + call compile_substr_fparser( this, & ! no compilation required + sub_str(2:sub_str_len-1) ) + return + + elseif (is_alpha(sub_str(1:1))) then + fidx = is_function( sub_str , f_len, n_params) + if (fidx > 0) then ! case 3: str = 'fcn(...)' + if (if_enclosed_parent(sub_str(f_len+1:))) then + i = f_len + 2 ! compile the parameters first + + do p = 2, functions(fidx)%n_params ! parse multiple parameters + j = comma_pos(sub_str(i:)) + + call compile_substr_fparser( this, & + sub_str(i:i+j-2) ) + i = i + j + enddo + + call compile_substr_fparser( this, & + sub_str(i:sub_str_len-1) ) + + ! add the compiled function to the code + call add_function( this, fidx ) + + ! decrease stack pointer if n_params > 1 + this%stack_ptr = this%stack_ptr - functions(fidx)%n_params + 1 + + return + endif + endif + + elseif (sub_str(1:1) == '-') then + if (if_enclosed_parent(sub_str(2:))) then ! case 4: str = '-(...)' + ! the same as neg(...) + call compile_substr_fparser( this, & ! compile the parameters first + sub_str(3:sub_str_len-1) ) + + ! add the compiled function to the code + call add_function( this, func_neg ) + + return + + elseif (is_alpha(sub_str(2:2))) then + + fidx = is_function( sub_str(2:) , f_len, n_params) + + if (fidx > 0) then ! case 5: str = '-fcn(...)' + if (if_enclosed_parent(sub_str(2+f_len:))) then + i = f_len + 3 ! compile the parameters first + + do p = 2, functions(fidx)%n_params ! parse multiple parameters + j = index(sub_str(i:),',') + call compile_substr_fparser( this, & + sub_str(i:i+j-2) ) + i = i + j + enddo + + call compile_substr_fparser( this, & + sub_str(i:sub_str_len-1) ) + + ! add the compiled function to the code followed by neg() + call add_function( this, fidx ) + call add_function( this, func_neg ) + + return + endif + endif + + endif + endif + + ! parse operators + ! only base level (k == 0) will be processed + ! the remaining levels are processed recursively + ! Operators are evaluated from left to right i.e. a - b - c ==> (a - b) - c + + ! Loop over operator groups + op_code = 0 + list_idx = 1 + do while ( list_idx <= size( oper ) ) + + n_op = same_priority_ops( list_idx ) + + k = 0 + + ! Process string from end to beginning + do i=sub_str_len,2,-1 + + ! skip sections inside parenthesis + select case (sub_str(i:i)) + case (')') + k = k+1 + case ('(') + k = k-1 + case default + continue + end select + + if (k == 0) then + + op_idx = scan_ops( sub_str, i, list_idx, n_op, op_len ) + + if ( op_idx > 0 ) then + + if ( num_signed_exp( i, sub_str ) ) then + ! number with signed exponent e.g. 1.0e-5 + ! this is a more complicated test because a variable name can end in + ! e or d + + ! ignore + continue + elseif ( (sub_str(1:1) == '-') .and. ( op_code >= op_div ) ) then + ! This is the only monadic operator we have + ! case 6: str = - ... op ... with op with higher priority than * or / (i.e. ^) + ! compile as neg(... op ...) + + ! compile the left operand + call compile_substr_fparser( this, sub_str(2:i-1) ) + + ! compile the right operand + call compile_substr_fparser( this, sub_str(i+op_len:) ) + + ! add the compiled operator to the code + call add_operator( this, op_code + op_idx ) + + ! decrease stack pointer + ! because the two values are converted into a single one + this%stack_ptr = this%stack_ptr - 1 + + call add_function( this, func_neg ) + + return + + else ! case 7: str = ... op ... + + ! compile the left operand + call compile_substr_fparser( this, sub_str(1:i-1) ) + + ! compile the right operand + call compile_substr_fparser( this, sub_str(i+op_len:) ) + + ! add the compiled operator to the code + call add_operator( this, op_code + op_idx ) + + ! decrease stack pointer + ! because the two values are converted into a single one + this%stack_ptr = this%stack_ptr - 1 + + return + endif + + endif + endif + + enddo + + ! Process next group of operators + list_idx = list_idx + n_op + 1 + op_code = op_code + n_op + enddo + + ! parse numbers and variables + + i = 1 + + ! check for var and -var + if (sub_str(1:1) == '-') then + varidx = is_variable_fparser( this, sub_str(2:), var_len) + if (varidx > 0) i = i + 1 + else + varidx = is_variable_fparser( this, sub_str(1:), var_len) + endif + + ! add the compiled value + if (varidx > 0) then ! variable + call add_variable( this, varidx ) + else ! number + call add_constant( this, strtodouble(sub_str, ierr) ) + endif + + ! increase stack pointer because one value will be added to the stack + this%stack_ptr = this%stack_ptr + 1 + + ! find maximum stack size + if (this%stack_ptr > this%stack_size) this%stack_size = this%stack_ptr + + ! if -var found add a neg() function + if (i>1) call add_function( this, func_neg ) + + contains + + !------------------------------------------------- + ! number with signed exponent e.g. 1.0e-5 + ! this is a more complicated test because a variable name can end in + ! e or d + !------------------------------------------------- + function num_signed_exp( i, sub_str ) + + implicit none + + integer, intent(in) :: i + character(len=*), intent(in) :: sub_str + logical :: num_signed_exp + + num_signed_exp = .false. + if ( i > 2 ) then + num_signed_exp = ((sub_str(i:i)=='-') .or. (sub_str(i:i)=='+')) .and. & + (scan( sub_str(i-1:i-1), 'edED' ) > 0) .and. & + (scan( sub_str(i-2:i-2), '0123456789.') > 0) + endif + end function num_signed_exp + !------------------------------------------------- + + !------------------------------------------------- + ! Check if str begins with any of the strings in oplist, and returns which one + !------------------------------------------------- + function scan_ops( in_str, i, list_idx, n_op, op_len ) + + implicit none + + character(len = *), intent(in) :: in_str + integer, intent(in) :: i, list_idx, n_op + integer, intent(out) :: op_len + + integer :: scan_ops, j, l + + do scan_ops = 1, n_op + j = list_idx + scan_ops - 1 + op_len = len_trim( oper( j ) ) + + l = min( op_len, len(in_str) - i + 1) + if (in_str(i:i+l-1) == oper(j)) return + enddo + + op_len = 0 + scan_ops = -1 + + end function scan_ops + !------------------------------------------------- + + !------------------------------------------------- + ! count operators with same priority + !------------------------------------------------- + function same_priority_ops( i ) + + implicit none + + integer, intent(in) :: i + integer :: same_priority_ops + + same_priority_ops = 0 + do while ( i + same_priority_ops <= size( oper ) ) + if ( oper(i + same_priority_ops) == ' ' ) exit + same_priority_ops = same_priority_ops + 1 + enddo + + end function same_priority_ops + !------------------------------------------------- + + +end subroutine compile_substr_fparser +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! add bytecode for a variable to the code stack +!----------------------------------------------------------------------------------------- + subroutine add_variable( this, varidx ) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + integer, intent(in) :: varidx + + ! local variables + + ! print *, "add variable = ", varidx + + this%code_size = this%code_size + 1 + this%bytecode(this%code_size) = bcode_var + varidx + + end subroutine add_variable +!----------------------------------------------------------------------------------------- + + +!----------------------------------------------------------------------------------------- +! add bytecode for a function to the code stack +!----------------------------------------------------------------------------------------- + subroutine add_function( this, fidx ) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + integer, intent(in) :: fidx + + ! local variables + + ! print *, "add func = '", trim(functions(fidx)%name),"'" + + this%code_size = this%code_size + 1 + this%bytecode(this%code_size) = n_oper + fidx + + end subroutine add_function +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! add bytecode for an operator to the code stack +!----------------------------------------------------------------------------------------- + subroutine add_operator( this, opidx ) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + integer, intent(in) :: opidx + + ! local variables + + ! print *, "add op = ",opidx + + this%code_size = this%code_size + 1 + this%bytecode(this%code_size) = opidx + + end subroutine add_operator +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! add bytecode for a number to the code stack +! and stores number in the intern_float stack +!----------------------------------------------------------------------------------------- + subroutine add_constant( this, float ) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + real(p_k_fparse), intent(in) :: float + + ! local variables + + !print *, 'add const = ', float + + this%code_size = this%code_size + 1 + this%data_size = this%data_size + 1 + this%data_stack(this%data_size) = float + this%bytecode(this%code_size) = bcode_data + + end subroutine add_constant +!----------------------------------------------------------------------------------------- + + +!----------------------------------------------------------------------------------------- +! test if current string is completely enclosed in parenthesis +!----------------------------------------------------------------------------------------- + function if_enclosed_parent(in_str) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + character (len = *), intent(in) :: in_str + logical :: if_enclosed_parent + + ! local variables + integer :: str_len + integer :: j,k + + ! executable statements + + str_len = len(in_str) + if_enclosed_parent = .false. + + if ((in_str(1:1) == '(') .and. (in_str(str_len:str_len) == ')')) then + k = 0 + do j = 2, str_len-1 + select case (in_str(j:j)) + case ('(') + k = k+1 + case (')') + k = k-1 + case default + continue + end select + if (k < 0) return + enddo + if (k == 0) if_enclosed_parent = .true. + endif + + end function if_enclosed_parent +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! returns the position of the first ',' at the first parenthesis level +!----------------------------------------------------------------------------------------- + function comma_pos(in_str) +!----------------------------------------------------------------------------------------- + implicit none + + ! dummy variables + + character (len = *), intent(in) :: in_str + integer :: comma_pos + + ! local variables + integer :: str_len, k + + ! executable statements + + str_len = len(in_str) + + k = 0 + do comma_pos = 1, str_len + select case (in_str(comma_pos:comma_pos)) + case ('(') + k = k+1 + case (')') + k = k-1 + case (',') + if (k == 0) return + case default + continue + end select + if (k < 0) return + enddo + + comma_pos = -1 + + end function comma_pos +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! check syntax on the supplied function +!----------------------------------------------------------------------------------------- + subroutine check_syntax_fparser( this, ierr ) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + type( t_fparser ), intent(inout) :: this + integer, intent (out) :: ierr + + !local variables + + integer, dimension(MAX_NESTED_FUNC) :: parent_count ! parethesis count at the current level + integer, dimension(MAX_NESTED_FUNC) :: param_count ! parameters count at the current level + integer, dimension(MAX_NESTED_FUNC) :: n_params ! number of parameters for the function + ! at the current level + + integer :: level ! function level being processed + + integer :: i, & ! current position in the string being parsed + fun_len, & ! length of the string being parsed + r_len, & ! length of the number processed + op_len, & ! length of the operator processed + f_len, & ! lengh of the function processed + var_len, & ! length of the variable processed + tmp_nparam + + character(len=1) :: c ! character being processed + + character(len=len_trim(this%func)) :: func ! local copy of the string to process + + integer :: numberType ! type of number being processed + + ! executable statements + + ! using this local variables ensures the length is correct + func = trim(this%func) + fun_len = len(func) + + ! iniatilize + parent_count = 0 + i = 1 + ierr = 0 + + level = 1 + +parse: do + if (i > fun_len) then + ierr = err_eof; call syntax_error(func,i,ierr); return + endif + + c = func(i:i) + ! check for leading + or - + if (c == '-' .or. c == '+') then + i = i+1 + if (i > fun_len) then + ierr = err_eof; return + endif + ! an operator must not follow + if (is_operator(func(i:)) > 0) then + ierr = err_mult_op; call syntax_error(func,i,ierr); return + endif + c = func(i:i) + endif + + ! check for functions + if (is_function(func(i:), f_len , tmp_nparam ) > 0) then + i = i + f_len + if (i > fun_len) then + ierr = err_eof; return + endif + c = func(i:i) + if ( c /= '(' ) then + ierr = err_fun_par; call syntax_error(func,i,ierr); return + endif + + i = i + 1 + level = level + 1 + parent_count(level) = 0 + param_count(level) = 0 + n_params(level) = tmp_nparam + cycle parse + endif + + ! check for opening parenthesis + if ( c == '(' ) then + ! if opening ( increase parent_count and go to the beginning of the loop + parent_count(level) = parent_count(level) + 1 + i = i+1 + cycle parse + endif + + + ! check for numbers and variables + if (scan(c, '0123456789.') > 0) then + ! constant + + call parsenumber( func(i:), numberType, r_len) + if (numberType < 0) then + ierr = err_inv_num; call syntax_error(func,i,ierr); return + endif + i = i + r_len + if (i > fun_len) exit + c = this%func(i:i) + elseif (is_variable_fparser(this, func(i:), var_len) > 0) then ! variable + ! variable + + i = i + var_len + if (i > fun_len) exit + c = this%func(i:i) + + endif + + ! check for closing parenthesis + do while (c == ')') + if ( i <= 1 ) then + ierr = err_parent; return + endif + + ! decrement parethesis count + parent_count(level) = parent_count(level) - 1 + + + ! check if closing too many ) + if (parent_count(level) < 0) then + ! check if closing function parenthesis + if (level > 1) then + !check if number of parameters is correct + if (param_count(level) /= n_params(level)-1) then + ierr = err_param; call syntax_error(func,i,ierr); return + endif + + level = level - 1 + i = i+ 1 + if (i > fun_len) exit + cycle parse + endif + ierr = err_parent;call syntax_error(func,i,ierr); return + endif + + if (this%func(i-1:i-1) == '(') then + ierr = err_parent; call syntax_error(func,i,ierr);return + endif + + i = i+1 + if (i > fun_len) exit + c = this%func(i:i) + enddo + + ! check for new parameter to function + if (c == ',') then + if ( i <= 1 ) then + ierr = err_syntax; call syntax_error(func,i,ierr); return + endif + + if ( (func(i-1:i-1) == '(') .or. (func(i-1:i-1) == ',') ) then + ierr = err_param; call syntax_error(func,i,ierr); return + endif + + if (parent_count(level) /= 0) then + ierr = err_parent; call syntax_error(func,i,ierr); return + endif + + param_count(level) = param_count(level) + 1 + if (param_count(level) >= n_params(level)) then + ierr = err_param; call syntax_error(func,i,ierr); return + endif + i = i + 1 + cycle parse + endif + + ! if we reach this point a valid operator must follow + if (i > fun_len) exit + if (is_operator(func(i:), op_len) > 0) then + i = i + op_len + if (i > fun_len) then + ierr = err_eof; call syntax_error(func,i,ierr); return + endif + ! an operator must not follow + if ( is_operator(func(i:)) > 0 ) then + ierr = err_mult_op; call syntax_error(func,i,ierr); return + endif + else + ierr = err_syntax; call syntax_error(func,i,ierr); return + endif + + enddo parse + + if ((parent_count(1) > 0) .or. (level > 1)) then + ierr = err_parent + call syntax_error(func,i,ierr); + endif + + end subroutine check_syntax_fparser +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! issue the appropriate error message +!----------------------------------------------------------------------------------------- + subroutine syntax_error(in_str, pos, ierr, nparam ) +!----------------------------------------------------------------------------------------- + + character (len=*) , intent (in) :: in_str + integer , intent(in) :: pos, ierr + integer , intent(in), optional :: nparam + +! character (len=pos-1) :: err_point + character (len = 64) :: err_msg + integer :: local_nparam + + if (present(nparam)) then + local_nparam = nparam + else + local_nparam = 1 + endif + +! err_point = " " + + write(0,'(A)') "Error compiling function:" + write(0,'(A)') in_str + write(0,'(A,A)') repeat('-',pos-1), '^' + + select case (ierr) + case ( err_eof ) + err_msg = "End of line reached" + case ( err_mult_op ) + err_msg = "Consecutive operators" + case ( err_inv_num ) + err_msg = "Invalid number format" + case ( err_parent ) + err_msg = "Parenthesis mismatch" + case ( err_syntax ) + err_msg = "Syntax error" + case ( err_fun_par ) + err_msg = "An opening parenthesis '(' must follow a function" + case ( err_param ) + write(err_msg,'(A,I0,A)') "Invalid number of parameters, ",local_nparam," expected" + case default + err_msg = "Syntax error" + end select + + write(0,'(A,A,I0)') trim(err_msg), ", col = ", pos + + + end subroutine syntax_error +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! checks if in_str begins with a variable +!----------------------------------------------------------------------------------------- + function is_variable_fparser( this, in_str , var_len) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + type( t_fparser ), intent(in) :: this + character ( len = * ), intent(in) :: in_str + integer :: is_variable_fparser + integer, optional, intent(out):: var_len + + ! local variables + integer :: i, l + + ! executable statements + + i = 1 + + do while (i <= this%nvars) + l = min(len_trim(this%var_list(i)), len(in_str)) + if (in_str(1:l) == trim(this%var_list(i))) then + is_variable_fparser = i + if (present(var_len)) var_len = len_trim(this%var_list(i)) + return + endif + i = i+1 + enddo + + is_variable_fparser = -1 + + end function is_variable_fparser +!----------------------------------------------------------------------------------------- + + +!----------------------------------------------------------------------------------------- +! checks if in_str begins with an operator and return the opcode +!----------------------------------------------------------------------------------------- + function is_operator( in_str , op_len) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + character ( len = * ), intent(in) :: in_str + integer :: is_operator + integer, optional, intent(out):: op_len + + ! local variables + integer :: i, l + + ! executable statements + + i = 1 + is_operator = 0 + + do while (i <= size( oper )) + if ( oper(i) /= ' ' ) then + is_operator = is_operator + 1 + l = min(len_trim(oper(i)), len(in_str)) + if (in_str(1:l) == oper(i)) then + if (present(op_len)) op_len = len_trim(oper(i)) + return + endif + endif + i = i+1 + enddo + + is_operator = -1 + + end function is_operator +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! checks if in_str begins with a function +!----------------------------------------------------------------------------------------- + function is_function( in_str , f_len, n_params) +!----------------------------------------------------------------------------------------- + + implicit none + + ! dummy variables + + character(len=*), intent(in) :: in_str + integer :: is_function + integer, optional,intent(out):: f_len + integer, optional,intent(out):: n_params + + ! local variables + integer :: i, l + + ! executable statements + + i = 1 + + do while (i <= n_functions) + l = min(len_trim(functions(i)%name), len(in_str)) + if (in_str(1:l) == functions(i)%name) then + is_function = i + if (present(f_len)) f_len = len_trim(functions(i)%name) + if (present(n_params)) n_params = functions(i)%n_params + return + endif + i = i+1 + enddo + + is_function = -1 + + end function is_function +!----------------------------------------------------------------------------------------- + +!----------------------------------------------------------------------------------------- +! return in_str with white spaces (spaces and tabs) removed +!----------------------------------------------------------------------------------------- + function remove_spaces( in_str ) +!----------------------------------------------------------------------------------------- + + character ( len = * ), intent(in) :: in_str + character (len = len_trim(in_str)) :: remove_spaces + + integer :: i, j + remove_spaces = "" + j = 1 + do i = 1, len(in_str) + ! if (in_str(i:i) .ne. " ") then + if (.not. is_blank(in_str(i:i))) then + remove_spaces(j:j) = in_str(i:i) + j = j +1 + endif + enddo + + end function remove_spaces +!----------------------------------------------------------------------------------------- + +!------------------------------------------------------------------------------- +subroutine parsenumber( str, numberType, numberLen ) +!------------------------------------------------------------------------------- +! parses the supplied string looking for a valid numeric value. If found +! returns the length of the string holding the number and the number type (1 - integer, +! 2 - float), otherwise returns length = 0, and type = -1 +!------------------------------------------------------------------------------- + + implicit none + + character(len=*), intent(in) :: str + integer, intent(out) :: numberLen, numberType + + character(len=len_trim(str)) :: value_text + character :: c, c2 + integer :: pos + logical :: has_decpoint, has_value, has_exponent, finished + + pos = 0 + has_decpoint = .false. + has_value = .false. + has_exponent = .false. + finished = .false. + value_text = '' + + ! check for signed number + c = str(1:1) + if ( c == '+' .or. c== '-' ) pos = pos+1 + + do + pos = pos + 1 + if ( pos > len( str ) ) then + finished = .true. + exit + endif + c = str(pos:pos) + + ! check for decimal point + if ( c == '.' ) then + if ( .not. has_decpoint ) then + has_decpoint = .true. + else + ! error, two decimal points found + finished = .true. + endif + + else if ( is_digit(c) ) then + has_value = .true. + + else if ( c=='e' .or. c=='E' .or. c=='d' .or. c=='D' ) then + if ( (.not. has_value) .or. has_exponent ) then + ! of exponent found before mantissa, or exponent symbol found inside exponent + finished = .true. + endif + has_exponent = .true. + has_decpoint = .true. ! decimal points are not allowed in the exponent + + ! check if next character is + or - + pos = pos+1 + if ( pos > len( str ) ) then + finished = .true. + exit + endif + c2 = str( pos:pos ) + if ((c2 == '+') .or. (c2 == '-')) then + value_text = trim(value_text)//c//c2 + cycle + else + pos = pos - 1 + endif + + has_value = .false. ! an integer must follow + else + ! invalid character found + ! finish processing + + finished = .true. + endif + + + if ( finished ) then + pos = pos-1 + exit + else + value_text = trim(value_text)//c + endif + + enddo + + if ( .not. has_value ) then + numberLen = 0 + numberType = -1 + return + endif + + numberLen = len_trim(value_text) + if ( has_exponent .or. has_decpoint ) then + numberType = 2 + else + numberType = 1 + endif + + +! print *, ">", trim(value_text), "< len = ", numberLen, ' type = ', numberType + +end subroutine parseNumber +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +function is_digit(c) +!------------------------------------------------------------------------------- +! test if character is a digit +!------------------------------------------------------------------------------- + implicit none + + character, intent(in) :: c + logical :: is_digit + +! if ((c >= '0') .and. (c <= '9')) then +! is_digit = .true. +! else +! is_digit = .false. +! endif + + is_digit = (c >= '0') .and. (c <= '9') + +end function is_digit +!------------------------------------------------------------------------------- + + !------------------------------------------------------------------------------- +function is_blank(c) +!------------------------------------------------------------------------------- +! test if character is blank (white space) +!------------------------------------------------------------------------------- + implicit none + + character, intent(in) :: c + logical :: is_blank + + if ((c == p_space) .or. (c == p_tab)) then + is_blank = .true. + else + is_blank = .false. + endif + +end function is_blank +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +function is_alpha(c) +!------------------------------------------------------------------------------- +! test if character is a letter +!------------------------------------------------------------------------------- + implicit none + + character, intent(in) :: c + logical :: is_alpha + +! if (((c >= 'a') .and. (c <= 'z')) .or. & +! ((c >= 'A') .and. (c <= 'Z'))) then +! is_alpha = .true. +! else +! is_alpha = .false. +! endif + + is_alpha = ((c >= 'a') .and. (c <= 'z')) .or. ((c >= 'A') .and. (c <= 'Z')) + +end function is_alpha +!------------------------------------------------------------------------------- + + +!------------------------------------------------------------------------------- +function strtodouble(str, ierr) +!------------------------------------------------------------------------------- +! convert string to integer +!------------------------------------------------------------------------------- + implicit none + + character(len=*), intent(in) :: str + integer, intent(out) :: ierr + real( kind(1.0d0) ) :: strtodouble + + read( str, *, iostat = ierr ) strtodouble + +end function strtodouble +!------------------------------------------------------------------------------- + + +end module m_fparser diff --git a/source/make.GF_OPENMPI b/source/make.GF_OPENMPI index b2f08f0..8a1ed43 100644 --- a/source/make.GF_OPENMPI +++ b/source/make.GF_OPENMPI @@ -1,11 +1,14 @@ # name of the compiler -JSON_LIB = /usr/local/jsonfortran-gnu-6.3.0/lib/ -JSON_INC = /usr/local/jsonfortran-gnu-6.3.0/lib/ +JSON_LIB = /home/tdali92/json-fortran-8.3.0/build/lib +JSON_INC = /home/tdali92/json-fortran-8.3.0/build/include + FC_GF_OPENMPI = mpif90 -c -I$(JSON_INC) -fopenmp -CC_GF_OPENMPI = mpicc -c -O3 -std=c99 +CC_GF_OPENMPI = gcc-12 -c -O3 -std=c99 -fopenmp LINKER_GF_OPENMPI = mpif90 -fopenmp -O3 +LDF = -L$(JSON_LIB) -ljsonfortran -lhdf5_fortran -lhdf5 -ldl -lz + OPTS_DBL_GF_OPENMPI = -O3 -fdefault-real-8 -fdefault-double-8 OPTS_SGL_GF_OPENMPI = -O3 @@ -13,9 +16,8 @@ FORMAT_FREE_GF_OPENMPI = -ffree-form FORMAT_FIXED_GF_OPENMPI = -ffixed-form # hdf5 libraries - -HDF5_LIB =/usr/local/hdf5/lib -HDF5_INC =/usr/local/hdf5/include +HDF5_LIB = /home/linuxbrew/.linuxbrew/Cellar/hdf5-mpi/1.12.2_1/lib +HDF5_INC = /home/linuxbrew/.linuxbrew/Cellar/hdf5-mpi/1.12.2_1/include HDF_LIBPATH_GF_OPENMPI= -L$(HDF5_LIB) -lhdf5_fortran -lhdf5_hl -lhdf5 -lz HDF_INCLUDE_PATH_GF_OPENMPI= -I$(HDF5_INC) -I${HDF5_LIB} -lhdf5hl_fortran -lhdf5_fortran -lhdf5_hl -lhdf5 @@ -28,6 +30,9 @@ MEM_OPTIONS_GF_OPENMPI = # options for linking # use static libs -STATIC_LINK_GF_OPENMPI = +STATIC_LINK_GF_OPENMPI = -L$(JSON_LIB) -ljsonfortran # use shared libs -SHARED_LINK_GF_OPENMPI = +SHARED_LINK_GF_OPENMPI = + + + diff --git a/source/makefile b/source/makefile index faf3c77..b425648 100644 --- a/source/makefile +++ b/source/makefile @@ -25,9 +25,10 @@ FORMAT_FIXED = ${FORMAT_FIXED_${SYS_SFX}} # Objects list OBJS_BASE = \ -parallel_class.o parallel_pipe_class.o perrors_class.o dtimer.o spect2d_class.o \ +parallel_class.o parallel_pipe_class.o perrors_class.o dtimer.o param.o spect2d_class.o \ spect3d_class.o \ hdf5io_class.o \ +fparser.o \ ufield2d_class.o ufield2d_lib77.o ufield2d_lib.o \ fft2d_lib.o fft2d_lib77.o fft2d_class.o \ fpois2d_lib.o fpois2d_lib77.o fpois2d_class.o \ @@ -80,7 +81,7 @@ beam3d_class.o : beam3d_class.f03 fdist3d_class.o field3d_class.o spect3d_class. species2d_class.o : species2d_class.f03 spect2d_class.o perrors_class.o \ parallel_pipe_class.o field2d_class.o hdf5io_class.o \ - fdist2d_class.o part2d_class.o + fdist2d_class.o part2d_class.o param.o ${FC} ${OPTS} ${FORMAT_FREE} species2d_class.f03 -o species2d_class.o @@ -165,6 +166,9 @@ ufield2d_class.o : ufield2d_class.f03 spect2d_class.o perrors_class.o \ hdf5io_class.o : hdf5io_class.f03 spect3d_class.o perrors_class.o parallel_pipe_class.o ${FC} ${OPTS} ${FORMAT_FREE} hdf5io_class.f03 -o hdf5io_class.o +fparser.o : fparser.f03 + ${FC} ${OPTS} ${FORMAT_FREE} fparser.f03 -o fparser.o + spect3d_class.o : spect3d_class.f03 perrors_class.o parallel_pipe_class.o spect2d_class.o ${FC} ${OPTS} ${FORMAT_FREE} spect3d_class.f03 -o spect3d_class.o @@ -174,6 +178,9 @@ spect2d_class.o : spect2d_class.f03 perrors_class.o parallel_pipe_class.o dtimer.o : dtimer.c ${CC} dtimer.c +param.o : param.f03 + ${FC} ${OPTS} ${FORMAT_FREE} param.f03 -o param.o + perrors_class.o : perrors_class.f03 parallel_class.o dtimer.o ${FC} ${OPTS} ${FORMAT_FREE} perrors_class.f03 -o perrors_class.o diff --git a/source/param.f03 b/source/param.f03 new file mode 100644 index 0000000..773960f --- /dev/null +++ b/source/param.f03 @@ -0,0 +1,15 @@ +module param + +implicit none + +public + + +! ================================================================ +! particle pusher +! ================================================================ +integer, parameter :: p_push2_std = 0, p_push2_robust = 1, p_push2_robust_subcyc = 2, & + p_push2_clamp = 3, p_push2_std_pgc = 4, p_push2_robust_pgc = 5 + + +end module param diff --git a/source/part2d_class.f03 b/source/part2d_class.f03 index bcbc286..b2afafd 100755 --- a/source/part2d_class.f03 +++ b/source/part2d_class.f03 @@ -28,6 +28,7 @@ module part2d_class use part2d_lib use hdf5io_class use mpi + use param implicit none @@ -74,6 +75,7 @@ module part2d_class generic :: del => end_part2d generic :: qdp => qdeposit generic :: amjdp => amjdeposit + generic :: amjdp_robust => amjdeposit_robust generic :: push => partpush generic :: pmv => pmove generic :: extpsi => extractpsi @@ -86,6 +88,7 @@ module part2d_class procedure, private :: end_part2d procedure, private :: qdeposit procedure, private :: amjdeposit + procedure, private :: amjdeposit_robust procedure, private :: partpush procedure, private :: pmove procedure, private :: extractpsi @@ -321,6 +324,38 @@ subroutine amjdeposit(this,ef,bf,psit,cu,amu,dcu,dex) call this%err%werrfl2(class//sname//' ended') end subroutine amjdeposit + + subroutine amjdeposit_robust(this,ef,bf,psit,cu,amu,dcu,dex) +! deposit the current, acceleration and momentum flux + + implicit none + + class(part2d), intent(inout) :: this + class(ufield2d), pointer, intent(in) :: cu, amu, dcu + class(ufield2d), pointer, intent(in) :: ef, bf, psit + real, intent(in) :: dex + character(len=18), save :: sname = 'amjdeposit_robust' +! local data + real, dimension(:,:,:), pointer :: pef => null(), pbf => null() + real, dimension(:,:,:), pointer :: ppsit => null(), pcu => null() + real, dimension(:,:,:), pointer :: pamu => null(), pdcu => null() + integer :: noff, nyp, nx, nxv, nypmx + + call this%err%werrfl2(class//sname//' started') + pef => ef%getrf(); pbf => bf%getrf() + ppsit => psit%getrf(); pcu => cu%getrf() + pamu => amu%getrf(); pdcu => dcu%getrf() + noff = ef%getnoff() + nxv = size(pef,2); nypmx = size(pef,3) + nx = ef%getnd1(); nyp = ef%getnd2p() + call PPGRDCJPPOST2L_QP_ROBUST(this%ppart,pef,pbf,ppsit(1,:,:),pcu,pdcu,& + &pamu,this%kpic,noff,nyp,this%qbm, this%dt,this%ci,this%xdim,& + &this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + + + call this%err%werrfl2(class//sname//' ended') + + end subroutine amjdeposit_robust ! subroutine partpush(this,ef,bf,psit,dex) @@ -608,13 +643,14 @@ subroutine pmove(this,fd) end subroutine pmove ! - subroutine extractpsi(this,psi,dex) + subroutine extractpsi(this,psi,dex, push_type) implicit none class(part2d), intent(inout) :: this class(ufield2d), pointer, intent(in) :: psi real, intent(in) :: dex + integer, intent(in) :: push_type character(len=18), save :: sname = 'extractpsi' ! local data real, dimension(:,:,:), pointer :: ppsi @@ -626,8 +662,15 @@ subroutine extractpsi(this,psi,dex) nyp = psi%getnd2p(); nx = psi%getnd1() nxv = size(ppsi,2); nypmx = size(ppsi,3) - call WPGPSIPOST2L_QP(this%ppart,ppsi(1,:,:),this%kpic,this%qbm,noff,& - &nyp,this%xdim,this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + select case (push_type) + case (p_push2_std) + call WPGPSIPOST2L_QP(this%ppart,ppsi(1,:,:),this%kpic,this%qbm,noff,& + &nyp,this%xdim,this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + case(p_push2_robust) + call WPGPSIPOST2L_QP_ROBUST(this%ppart,ppsi(1,:,:),this%kpic,this%qbm,noff,& + &nyp,this%xdim,this%nppmx0,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + end select + call this%err%werrfl2(class//sname//' ended') diff --git a/source/part2d_lib.f03 b/source/part2d_lib.f03 index 1ac04cf..bfb46b8 100644 --- a/source/part2d_lib.f03 +++ b/source/part2d_lib.f03 @@ -89,6 +89,26 @@ subroutine PPGRDCJPPOST2L_QP(ppart,fxy,bxy,psit,cu,dcu,amu,kpic,noff,& real, intent(in) :: dex end subroutine end interface + ! + interface + subroutine PPGRDCJPPOST2L_QP_ROBUST(ppart,fxy,bxy,psit,cu,dcu,amu,kpic,noff,& + &nyp,qbm,dt,ci,idimp,nppmx,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + implicit none + integer, intent(in) :: noff, nyp, idimp, nppmx, nx, mx, my + integer, intent(in) :: nxv, nypmx, mx1, mxyp1 + real, intent(in) :: qbm, dt, ci + real, dimension(idimp,nppmx,mxyp1), intent(inout) :: ppart + real, dimension(2,nxv,nypmx), intent(in) :: fxy + real, dimension(3,nxv,nypmx), intent(in) :: bxy + real, dimension(nxv,nypmx), intent(in) :: psit + real, dimension(3,nxv,nypmx), intent(inout) :: cu + real, dimension(2,nxv,nypmx), intent(inout) :: dcu + real, dimension(3,nxv,nypmx), intent(inout) :: amu + integer, dimension(mxyp1), intent(in) :: kpic + real, intent(in) :: dex + end subroutine + end interface + ! interface subroutine PPGRBPPUSH23L_QP(ppart,fxy,bxy,psit,kpic,noff,nyp, & @@ -207,6 +227,19 @@ subroutine PPPRSTOR2L(ppart,ppbuff,ncl,ihole,idimp,nppmx,mxyp1,& integer, dimension(2,ntmax+1,mxyp1), intent(in) :: ihole end subroutine end interface +! + interface + subroutine WPGPSIPOST2L_QP_ROBUST(ppart,psi,kpic,qbm,noff,nyp,idimp,np& + &pmx,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + implicit none + integer, intent(in) :: noff, nyp, idimp, nppmx, nx, mx, my + integer, intent(in) :: mx1, mxyp1, nxv, nypmx + real, intent(in) :: dex,qbm + real, dimension(idimp,nppmx,mxyp1), intent(inout) :: ppart + real, dimension(nxv,nypmx), intent(in) :: psi + integer, dimension(mxyp1), intent(in) :: kpic + end subroutine + end interface ! interface subroutine WPGPSIPOST2L_QP(ppart,psi,kpic,qbm,noff,nyp,idimp,np& diff --git a/source/part2d_lib77.f b/source/part2d_lib77.f index 2d6ff8d..35a39c2 100644 --- a/source/part2d_lib77.f +++ b/source/part2d_lib77.f @@ -768,6 +768,371 @@ subroutine PPGRDCJPPOST2L_QP(ppart,fxy,bxy,psit,cu,dcu,amu,kpic,no !$OMP END PARALLEL DO return end + +!----------------------------------------------------------------------- + subroutine PPGRDCJPPOST2L_QP_ROBUST(ppart,fxy,bxy,psit,cu,dcu,amu,& + &kpic,noff,nyp,qbm,dt,ci,idimp,nppmx,nx,mx,my,nxv,nypmx,mx1,mxyp1,d& + &ex) +! for 2-1/2d code, this subroutine calculates particle momentum flux, +! acceleration density and current density using first-order spline +! interpolation for relativistic particles. + implicit none + integer noff, nyp, idimp, nppmx, nx, mx, my, nxv, nypmx + integer mx1, mxyp1 + real qbm, dt, ci,dex + real ppart, fxy, bxy, cu, dcu, amu, psit + integer kpic + dimension ppart(idimp,nppmx,mxyp1) + dimension fxy(2,nxv,nypmx), bxy(3,nxv,nypmx), psit(nxv,nypmx) + dimension cu(3,nxv,nypmx), dcu(2,nxv,nypmx), amu(3,nxv,nypmx) + dimension kpic(mxyp1) +! local data +! integer MXV, MYV +! parameter(MXV=33,MYV=33) + integer noffp, moffp, nppp + integer mnoff, i, j, k, nn, mm + real qtmh, dti, ci2, gami, qtmg, gh, dxp, dyp, amx, amy + real dx, dy, dz, ox, oy, oz, qm + real acx, acy, acz, omxt, omyt, omzt, omt, anorm + real rot1, rot2, rot3, rot4, rot5, rot6, rot7, rot8, rot9 + real x, y, vx, vy, vz, p2, v1, v2, v3, v4 + real sfxy, sbxy, scu, sdcu, samu, spsit + real tmpx,tmpy,tmpz + real dotprod + dimension sfxy(2,mx+1,my+1), sbxy(3,mx+1,my+1) + dimension spsit(mx+1,my+1) + dimension scu(3,mx+1,my+1), sdcu(2,mx+1,my+1), samu(3,mx+1,my+1) +! dimension sfxy(2,MXV,MYV), sbxy(3,MXV,MYV) +! dimension spsit(MXV,MYV) +! dimension scu(3,MXV,MYV), sdcu(2,MXV,MYV), samu(3,MXV,MYV) + real qm1, qtmh1,qtmh2,idex,inv_part_7,ddx,ddy,ddz,p6,p7 +! dimension sfxy(3,mx+1,my+1), sbxy(3,mx+1,my+1) +! dimension scu(3,mx+1,my+1,), sdcu(3,mx+1,my+1), samu(4,mx+1,my+1) + qtmh = 0.5*qbm*dt + dti = 1.0/dt + ci2 = ci*ci + idex = 1.0/dex +! error if local array is too small +! if ((mx.ge.MXV).or.(my.ge.MYV)) return +! loop over tiles +!$OMP PARALLEL DO +!$OMP& PRIVATE(i,j,k,noffp,moffp,nppp,nn,mm,mnoff,x,y,vx,vy,vz, +!$OMP& v1,v2,v3,dxp,dyp,amx,amy,dx,dy,dz,ox,oy,oz,acx,acy,acz, +!$OMP& omxt,omyt,omzt,omt,anorm,rot1,rot2,rot3,rot4,rot5, +!$OMP& tmpx,tmpy,tmpz, +!$OMP& rot6,rot7,rot8,rot9,spsit,sfxy,sbxy,scu,sdcu,samu,qm1, +!$OMP& qtmh1,qtmh2,inv_part_7,ddx,ddy,ddz,p6,p7,qm, dotprod) + do 120 k = 1, mxyp1 + noffp = (k - 1)/mx1 + moffp = my*noffp + noffp = mx*(k - mx1*noffp - 1) + nppp = kpic(k) + mnoff = moffp + noff - 1 +! load local fields from global arrays + nn = min(mx,nx-noffp) + 1 + mm = min(my,nyp-moffp) + 1 + do 20 j = 1, mm + do 10 i = 1, nn + sfxy(1,i,j) = fxy(1,i+noffp,j+moffp) + sfxy(2,i,j) = fxy(2,i+noffp,j+moffp) + 10 continue + 20 continue + do 40 j = 1, mm + do 30 i = 1, nn + sbxy(1,i,j) = bxy(1,i+noffp,j+moffp) + sbxy(2,i,j) = bxy(2,i+noffp,j+moffp) + sbxy(3,i,j) = bxy(3,i+noffp,j+moffp) + ! print *, i+noffp, j+moffp, sbxy(1:3,i,j) + 30 continue + 40 continue + do 46 j = 1, mm + do 42 i = 1, nn + spsit(i,j) = psit(i+noffp,j+moffp) + 42 continue + 46 continue +! zero out local accumulators + do 60 j = 1, my+1 + do 50 i = 1, mx+1 + samu(1,i,j) = 0.0 + samu(2,i,j) = 0.0 + samu(3,i,j) = 0.0 + sdcu(1,i,j) = 0.0 + sdcu(2,i,j) = 0.0 + scu(1,i,j) = 0.0 + scu(2,i,j) = 0.0 + scu(3,i,j) = 0.0 + 50 continue + 60 continue +! loop over particles in tile + do 70 j = 1, nppp +! find interpolation weights + x = ppart(1,j,k) + y = ppart(2,j,k) + p6 = ppart(6,j,k) + p7 = ppart(7,j,k) + qm = ppart(8,j,k) + nn = x + mm = y + dxp = x - real(nn) + dyp = y - real(mm) + nn = nn - noffp + 1 + mm = mm - mnoff + amx = 1.0 - dxp + amy = 1.0 - dyp + + inv_part_7 = 1.0 / p7 + qtmh1 = qtmh * inv_part_7 + qtmh2 = qtmh1 * p6 + +! find electric field + dx = amx*sfxy(1,nn,mm) + dy = amx*sfxy(2,nn,mm) + dz = amx*spsit(nn,mm) + + dx = amy*(dxp*sfxy(1,nn+1,mm) + dx) + dy = amy*(dxp*sfxy(2,nn+1,mm) + dy) + dz = amy*(dxp*spsit(nn+1,mm) + dz) + + acx = amx*sfxy(1,nn,mm+1) + acy = amx*sfxy(2,nn,mm+1) + acz = amx*spsit(nn,mm+1) + + dx = dx + dyp*(dxp*sfxy(1,nn+1,mm+1) + acx) + dy = dy + dyp*(dxp*sfxy(2,nn+1,mm+1) + acy) + dz = dz + dyp*(dxp*spsit(nn+1,mm+1) + acz) + +! find magnetic field + ox = amx*sbxy(1,nn,mm) + oy = amx*sbxy(2,nn,mm) + oz = amx*sbxy(3,nn,mm) + + ox = amy*(dxp*sbxy(1,nn+1,mm) + ox) + oy = amy*(dxp*sbxy(2,nn+1,mm) + oy) + oz = amy*(dxp*sbxy(3,nn+1,mm) + oz) + + acx = amx*sbxy(1,nn,mm+1) + acy = amx*sbxy(2,nn,mm+1) + acz = amx*sbxy(3,nn,mm+1) + + ox = ox + dyp*(dxp*sbxy(1,nn+1,mm+1) + acx) + oy = oy + dyp*(dxp*sbxy(2,nn+1,mm+1) + acy) + oz = oz + dyp*(dxp*sbxy(3,nn+1,mm+1) + acz) + +!calculate half impulse + ddx = qtmh2*(oy - dx) + ddy = qtmh2*(- ox - dy) + ddz = qtmh2* dz * idex + + ox = ox * dex + oy = oy * dex + + vx = ppart(3,j,k) + vy = ppart(4,j,k) + vz = ppart(5,j,k) + +! half acceleration + acx = vx + ddx + acy = vy + ddy + acz = vz + ddz + +! calculate cyclotron frequency + omxt = qtmh1*ox + omyt = qtmh1*oy + omzt = qtmh1*oz +! calculate rotation matrix + dotprod = (acx * omxt + acy * omyt + acz * omzt) + omt = omxt**2 + omyt**2 + omzt**2 + anorm = 2.0/(1.0 + omt) + omt = 0.5*(1.0 - omt) +! new momentum + v1 = (omt * acx+ dotprod * omxt + acy*omzt-acz*omyt)*anorm + ddx + v2 = (omt * acy+ dotprod * omyt + acz*omxt-acx*omzt)*anorm + ddy + v3 = (omt * acz+ dotprod * omzt + acx*omyt-acy*omxt)*anorm + ddz + +! deposit momentum flux, acceleration density, and current density + + ox = 0.5*(v1 + vx) + oy = 0.5*(v2 + vy) + oz = 0.5*(v3 + vz) + ppart(6,j,k) = sqrt(1.0 + (ox* ox + oy* oy + oz* oz) * dex * dex) + ppart(7,j,k) = ppart(6,j,k) - oz * dex + inv_part_7 = 1.0/ppart(7,j,k) + + qm1 = qm*inv_part_7 + + amx = qm1*amx + dxp = qm1*dxp + + + vx = (v1 - vx)*dti + vy = (v2 - vy)*dti + + dz = qbm*(dz + (dx*ox+dy*oy)*dex*dex*inv_part_7) + + + dx = amx*amy + dy = dxp*amy + + v1 = ox*ox*inv_part_7 + v2 = oy*oy*inv_part_7 + v3 = ox*oy*inv_part_7 + + vx = vx+ox*dz*inv_part_7 + vy = vy+oy*dz*inv_part_7 + + samu(1,nn,mm) = samu(1,nn,mm) + v1*dx + samu(2,nn,mm) = samu(2,nn,mm) + v2*dx + samu(3,nn,mm) = samu(3,nn,mm) + v3*dx + + sdcu(1,nn,mm) = sdcu(1,nn,mm) + vx*dx + sdcu(2,nn,mm) = sdcu(2,nn,mm) + vy*dx + + scu(1,nn,mm) = scu(1,nn,mm) + ox*dx + scu(2,nn,mm) = scu(2,nn,mm) + oy*dx + scu(3,nn,mm) = scu(3,nn,mm) + oz*dx + + dx = amx*dyp + samu(1,nn+1,mm) = samu(1,nn+1,mm) + v1*dy + samu(2,nn+1,mm) = samu(2,nn+1,mm) + v2*dy + samu(3,nn+1,mm) = samu(3,nn+1,mm) + v3*dy + + sdcu(1,nn+1,mm) = sdcu(1,nn+1,mm) + vx*dy + sdcu(2,nn+1,mm) = sdcu(2,nn+1,mm) + vy*dy + + scu(1,nn+1,mm) = scu(1,nn+1,mm) + ox*dy + scu(2,nn+1,mm) = scu(2,nn+1,mm) + oy*dy + scu(3,nn+1,mm) = scu(3,nn+1,mm) + oz*dy + + dy = dxp*dyp + samu(1,nn,mm+1) = samu(1,nn,mm+1) + v1*dx + samu(2,nn,mm+1) = samu(2,nn,mm+1) + v2*dx + samu(3,nn,mm+1) = samu(3,nn,mm+1) + v3*dx + + sdcu(1,nn,mm+1) = sdcu(1,nn,mm+1) + vx*dx + sdcu(2,nn,mm+1) = sdcu(2,nn,mm+1) + vy*dx + + scu(1,nn,mm+1) = scu(1,nn,mm+1) + ox*dx + scu(2,nn,mm+1) = scu(2,nn,mm+1) + oy*dx + scu(3,nn,mm+1) = scu(3,nn,mm+1) + oz*dx + + samu(1,nn+1,mm+1) = samu(1,nn+1,mm+1) + v1*dy + samu(2,nn+1,mm+1) = samu(2,nn+1,mm+1) + v2*dy + samu(3,nn+1,mm+1) = samu(3,nn+1,mm+1) + v3*dy + + sdcu(1,nn+1,mm+1) = sdcu(1,nn+1,mm+1) + vx*dy + sdcu(2,nn+1,mm+1) = sdcu(2,nn+1,mm+1) + vy*dy + + scu(1,nn+1,mm+1) = scu(1,nn+1,mm+1) + ox*dy + scu(2,nn+1,mm+1) = scu(2,nn+1,mm+1) + oy*dy + scu(3,nn+1,mm+1) = scu(3,nn+1,mm+1) + oz*dy + 70 continue +! deposit currents to interior points in global array + nn = min(mx,nxv-noffp) + mm = min(my,nypmx-moffp) + do 90 j = 2, mm + do 80 i = 2, nn + amu(1,i+noffp,j+moffp) = amu(1,i+noffp,j+moffp) + samu(1,i,j) + amu(2,i+noffp,j+moffp) = amu(2,i+noffp,j+moffp) + samu(2,i,j) + amu(3,i+noffp,j+moffp) = amu(3,i+noffp,j+moffp) + samu(3,i,j) + + dcu(1,i+noffp,j+moffp) = dcu(1,i+noffp,j+moffp) + sdcu(1,i,j) + dcu(2,i+noffp,j+moffp) = dcu(2,i+noffp,j+moffp) + sdcu(2,i,j) + + cu(1,i+noffp,j+moffp) = cu(1,i+noffp,j+moffp) + scu(1,i,j) + cu(2,i+noffp,j+moffp) = cu(2,i+noffp,j+moffp) + scu(2,i,j) + cu(3,i+noffp,j+moffp) = cu(3,i+noffp,j+moffp) + scu(3,i,j) + 80 continue + 90 continue +! deposit currents to edge points in global array + mm = min(my+1,nypmx-moffp) + do 100 i = 2, nn +!$OMP ATOMIC + amu(1,i+noffp,1+moffp) = amu(1,i+noffp,1+moffp) + samu(1,i,1) +!$OMP ATOMIC + amu(2,i+noffp,1+moffp) = amu(2,i+noffp,1+moffp) + samu(2,i,1) +!$OMP ATOMIC + amu(3,i+noffp,1+moffp) = amu(3,i+noffp,1+moffp) + samu(3,i,1) +!$OMP ATOMIC + dcu(1,i+noffp,1+moffp) = dcu(1,i+noffp,1+moffp) + sdcu(1,i,1) +!$OMP ATOMIC + dcu(2,i+noffp,1+moffp) = dcu(2,i+noffp,1+moffp) + sdcu(2,i,1) +!$OMP ATOMIC + cu(1,i+noffp,1+moffp) = cu(1,i+noffp,1+moffp) + scu(1,i,1) +!$OMP ATOMIC + cu(2,i+noffp,1+moffp) = cu(2,i+noffp,1+moffp) + scu(2,i,1) +!$OMP ATOMIC + cu(3,i+noffp,1+moffp) = cu(3,i+noffp,1+moffp) + scu(3,i,1) + if (mm > my) then +!$OMP ATOMIC + amu(1,i+noffp,mm+moffp) = amu(1,i+noffp,mm+moffp) & + & + samu(1,i,mm) +!$OMP ATOMIC + amu(2,i+noffp,mm+moffp) = amu(2,i+noffp,mm+moffp) & + & + samu(2,i,mm) +!$OMP ATOMIC + amu(3,i+noffp,mm+moffp) = amu(3,i+noffp,mm+moffp) & + & + samu(3,i,mm) +!$OMP ATOMIC + dcu(1,i+noffp,mm+moffp) = dcu(1,i+noffp,mm+moffp) & + & + sdcu(1,i,mm) +!$OMP ATOMIC + dcu(2,i+noffp,mm+moffp) = dcu(2,i+noffp,mm+moffp) & + & + sdcu(2,i,mm) +!$OMP ATOMIC + cu(1,i+noffp,mm+moffp) = cu(1,i+noffp,mm+moffp) + scu(1,i,mm) +!$OMP ATOMIC + cu(2,i+noffp,mm+moffp) = cu(2,i+noffp,mm+moffp) + scu(2,i,mm) +!$OMP ATOMIC + cu(3,i+noffp,mm+moffp) = cu(3,i+noffp,mm+moffp) + scu(3,i,mm) + endif + 100 continue + nn = min(mx+1,nxv-noffp) + do 110 j = 1, mm +!$OMP ATOMIC + amu(1,1+noffp,j+moffp) = amu(1,1+noffp,j+moffp) + samu(1,1,j) +!$OMP ATOMIC + amu(2,1+noffp,j+moffp) = amu(2,1+noffp,j+moffp) + samu(2,1,j) +!$OMP ATOMIC + amu(3,1+noffp,j+moffp) = amu(3,1+noffp,j+moffp) + samu(3,1,j) +!$OMP ATOMIC + dcu(1,1+noffp,j+moffp) = dcu(1,1+noffp,j+moffp) + sdcu(1,1,j) +!$OMP ATOMIC + dcu(2,1+noffp,j+moffp) = dcu(2,1+noffp,j+moffp) + sdcu(2,1,j) +!$OMP ATOMIC + cu(1,1+noffp,j+moffp) = cu(1,1+noffp,j+moffp) + scu(1,1,j) +!$OMP ATOMIC + cu(2,1+noffp,j+moffp) = cu(2,1+noffp,j+moffp) + scu(2,1,j) +!$OMP ATOMIC + cu(3,1+noffp,j+moffp) = cu(3,1+noffp,j+moffp) + scu(3,1,j) + if (nn > mx) then +!$OMP ATOMIC + amu(1,nn+noffp,j+moffp) = amu(1,nn+noffp,j+moffp) & + & + samu(1,nn,j) +!$OMP ATOMIC + amu(2,nn+noffp,j+moffp) = amu(2,nn+noffp,j+moffp) & + & + samu(2,nn,j) +!$OMP ATOMIC + amu(3,nn+noffp,j+moffp) = amu(3,nn+noffp,j+moffp) & + & + samu(3,nn,j) +!$OMP ATOMIC + dcu(1,nn+noffp,j+moffp) = dcu(1,nn+noffp,j+moffp) & + & + sdcu(1,nn,j) +!$OMP ATOMIC + dcu(2,nn+noffp,j+moffp) = dcu(2,nn+noffp,j+moffp) & + & + sdcu(2,nn,j) +!$OMP ATOMIC + cu(1,nn+noffp,j+moffp) = cu(1,nn+noffp,j+moffp) + scu(1,nn,j) +!$OMP ATOMIC + cu(2,nn+noffp,j+moffp) = cu(2,nn+noffp,j+moffp) + scu(2,nn,j) +!$OMP ATOMIC + cu(3,nn+noffp,j+moffp) = cu(3,nn+noffp,j+moffp) + scu(3,nn,j) + endif + 110 continue + 120 continue +!$OMP END PARALLEL DO + return + end +!----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine PPGRBPPUSH23L_QP(ppart,fxy,bxy,psit,kpic,noff,nyp,qbm, & &dt,ci,ek,idimp,nppmx,nx,ny,mx,my,nxv,nypmx,mx1,mxyp1,ipbc,dex) @@ -2204,6 +2569,89 @@ subroutine PPPRSTOR2L(ppart,ppbuff,ncl,ihole,idimp,nppmx,mxyp1, & !$OMP END PARALLEL DO return end + +!----------------------------------------------------------------------- +!-deposit psi on particles + subroutine WPGPSIPOST2L_QP_ROBUST(ppart,psi,kpic,qbm,noff,nyp + 1,idimp,nppmx,nx,mx,my,nxv,nypmx,mx1,mxyp1,dex) + + implicit none + integer noff, nyp, idimp, nppmx, nx, mx, my, nxv, nypmx + integer mx1, mxyp1 + real dex,qbm + real ppart, psi + integer kpic + dimension ppart(idimp,nppmx,mxyp1) + dimension psi(nxv,nypmx) + dimension kpic(mxyp1) +! local data +! integer MXV, MYV +! parameter(MXV=33,MYV=33) + integer noffp, moffp, nppp + integer mnoff, i, j, k, nn, mm + real dxp, dyp, amx, amy, acx + real dx, dy, dz + real x, y, vx, vy, vz + real spsi + dimension spsi(mx+1,my+1) + real dx2 +! dimension sfxy(3,mx+1,my+1), sbxy(3,mx+1,my+1) +! dimension scu(3,mx+1,my+1,), sdcu(3,mx+1,my+1), samu(4,mx+1,my+1) + dx2 = dex * dex +! error if local array is too small +! if ((mx.ge.MXV).or.(my.ge.MYV)) return +! loop over tiles +!$OMP PARALLEL DO +!$OMP& PRIVATE(i,j,k,noffp,moffp,nppp,nn,mm,mnoff,x,y,vx,vy,vz, +!$OMP& dxp,dyp,amx,amy,dx,acx,spsi) + do 120 k = 1, mxyp1 + noffp = (k - 1)/mx1 + moffp = my*noffp + noffp = mx*(k - mx1*noffp - 1) + nppp = kpic(k) + mnoff = moffp + noff - 1 +! load local fields from global arrays + nn = min(mx,nx-noffp) + 1 + mm = min(my,nyp-moffp) + 1 + do 20 j = 1, mm + do 10 i = 1, nn + spsi(i,j) = psi(i+noffp,j+moffp) + 10 continue + 20 continue +! loop over particles in tile + do 70 j = 1, nppp +! find interpolation weights + x = ppart(1,j,k) + y = ppart(2,j,k) + vx = ppart(3,j,k) + vy = ppart(4,j,k) + vz = ppart(5,j,k) + nn = x + mm = y + dxp = x - real(nn) + dyp = y - real(mm) + nn = nn - noffp + 1 + mm = mm - mnoff + amx = 1.0 - dxp + amy = 1.0 - dyp + + +! find electric field + dx = amx*spsi(nn,mm) + dx = amy*(dxp*spsi(nn+1,mm) + dx) + acx = amx*spsi(nn,mm+1) + dx = dx + dyp*(dxp*spsi(nn+1,mm+1) + acx) + + dx = - dx*dx2*qbm + ppart(6,j,k) = sqrt(1.0 + dx2 * (vx* vx + vy* vy + vz* vz)) + ppart(7,j,k) = ppart(6,j,k) - dex * vz + + 70 continue + 120 continue +!$OMP END PARALLEL DO + return + end +!----------------------------------------------------------------------- !----------------------------------------------------------------------- !-deposit psi on particles subroutine WPGPSIPOST2L_QP(ppart,psi,kpic,qbm,noff,nyp,idimp,nppmx diff --git a/source/simulation_class.f03 b/source/simulation_class.f03 index 212ac9e..047ffaa 100644 --- a/source/simulation_class.f03 +++ b/source/simulation_class.f03 @@ -17,6 +17,7 @@ module simulation_class use hdf5io_class use input_class use mpi + use param implicit none @@ -438,6 +439,9 @@ subroutine init_sim_beams(this,input) case (3) allocate(fdist3d_003::this%pf(i)%p) call this%pf(i)%p%new(input,i) + case (13) + allocate(fdist3d_013::this%pf(i)%p) + call this%pf(i)%p%new(input,i) case (100) allocate(fdist3d_100::this%pf(i)%p) call this%pf(i)%p%new(input,i) @@ -511,9 +515,10 @@ subroutine init_sim_species(this,input,s) integer :: npf real, dimension(3,100) :: arg character(len=20) :: sn,s1 - integer :: indz, xppc, yppc + integer :: indz, xppc, yppc, push_type real :: min, max, cwp, n0 real :: qm, qbm, dz + character(len=:), allocatable :: str this%err => input%err this%p => input%pp @@ -539,6 +544,23 @@ subroutine init_sim_species(this,input,s) write (sn,'(I3.3)') i s1 = 'species('//trim(sn)//')' call input%get(trim(s1)//'.profile',npf) + + push_type = p_push2_std + + if (input%found(trim(s1)//'.push_type')) then + call input%get( trim(s1)//'.push_type', str ) + select case (trim(str)) + case ( 'standard' ) + push_type = p_push2_std + case ( 'robust' ) + push_type = p_push2_robust + case default + write (erstr,*) 'Invalid pusher type! Only "standard", "robust" & + &are supported currently.' + end select + + endif + select case (npf) case (0) allocate(fdist2d_000::this%pf(i)%p) @@ -552,6 +574,9 @@ subroutine init_sim_species(this,input,s) case (12) allocate(fdist2d_012::this%pf(i)%p) call this%pf(i)%p%new(input,i) + case (13) + allocate(fdist2d_013::this%pf(i)%p) + call this%pf(i)%p%new(input,i) ! Add new distributions right above this line case default write (erstr,*) 'Invalid species profile number:', npf @@ -562,7 +587,7 @@ subroutine init_sim_species(this,input,s) 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) + &qbm=qbm,dt=dz,ci=1.0,xdim=8,s=s,push_type=push_type) end do call this%err%werrfl2(class//sname//' ended') diff --git a/source/species2d_class.f03 b/source/species2d_class.f03 index 3fbd04b..f5d32ed 100644 --- a/source/species2d_class.f03 +++ b/source/species2d_class.f03 @@ -12,7 +12,7 @@ module species2d_class use field3d_class use part2d_class use hdf5io_class - + use param implicit none private @@ -31,6 +31,7 @@ module species2d_class class(field2d), pointer :: amu => null(), dcu => null() class(field3d), pointer :: q3 => null() class(fdist2d), pointer :: pf => null() + integer :: push_type contains @@ -71,7 +72,7 @@ module species2d_class contains ! - subroutine init_species2d(this,pp,perr,psp,pf,qbm,dt,ci,xdim,s) + subroutine init_species2d(this,pp,perr,psp,pf,qbm,dt,ci,xdim,s,push_type) implicit none @@ -82,6 +83,7 @@ subroutine init_species2d(this,pp,perr,psp,pf,qbm,dt,ci,xdim,s) class(fdist2d), intent(inout), target :: pf real, intent(in) :: qbm, dt, ci, s integer, intent(in) :: xdim + integer, intent(in) :: push_type ! local data character(len=18), save :: sname = 'init_species2d:' @@ -90,6 +92,7 @@ subroutine init_species2d(this,pp,perr,psp,pf,qbm,dt,ci,xdim,s) this%err => perr this%p => pp this%pf => pf + this%push_type = push_type call this%err%werrfl2(class//sname//' started') @@ -201,8 +204,14 @@ subroutine amjdp_species2d(this,ef,bf,psit,cu,amu,dcu,dex) call this%dcu%as(0.0) call this%amu%as(0.0) - call this%pd%amjdp(ef%getrs(),bf%getrs(),psit%getrs(),this%cu%getrs(),& - &this%amu%getrs(),this%dcu%getrs(),dex) + select case (this%push_type) + case (p_push2_std) + call this%pd%amjdp(ef%getrs(),bf%getrs(),psit%getrs(),this%cu%getrs(),& + &this%amu%getrs(),this%dcu%getrs(),dex) + case ( p_push2_robust ) + call this%pd%amjdp_robust(ef%getrs(),bf%getrs(),psit%getrs(),this%cu%getrs(),& + &this%amu%getrs(),this%dcu%getrs(),dex) + end select call this%cu%ag() call this%dcu%ag() @@ -265,7 +274,7 @@ subroutine extpsi_species2d(this,psi,dex) call this%err%werrfl2(class//sname//' started') - call this%pd%extpsi(psi%getrs(),dex) + call this%pd%extpsi(psi%getrs(),dex, this%push_type) call this%err%werrfl2(class//sname//' ended')