diff options
Diffstat (limited to 'Source/Initialization')
| -rw-r--r-- | Source/Initialization/Make.package | 1 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.H | 9 | ||||
| -rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 36 | ||||
| -rw-r--r-- | Source/Initialization/WarpX_parser.F90 | 1365 |
4 files changed, 17 insertions, 1394 deletions
diff --git a/Source/Initialization/Make.package b/Source/Initialization/Make.package index 829b9a91c..b825924c6 100644 --- a/Source/Initialization/Make.package +++ b/Source/Initialization/Make.package @@ -3,7 +3,6 @@ CEXE_sources += WarpXInitData.cpp CEXE_sources += CustomMomentumProb.cpp CEXE_sources += PlasmaInjector.cpp CEXE_headers += PlasmaInjector.H -F90EXE_sources += WarpX_parser.F90 INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Initialization VPATH_LOCATIONS += $(WARPX_HOME)/Source/Initialization diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index a25d1141b..7e61fa5d2 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -6,6 +6,7 @@ #include "AMReX_REAL.H" #include <AMReX_Vector.H> #include <WarpXConst.H> +#include <WarpXParser.H> #include "AMReX_ParmParse.H" #include "AMReX_Utility.H" @@ -70,7 +71,7 @@ public: UserConstants my_constants; private: std::string _parse_density_function; - int parser_instance_number=0; + WarpXParser parser_density; }; /// @@ -177,9 +178,9 @@ private: std::string _parse_momentum_function_ux; std::string _parse_momentum_function_uy; std::string _parse_momentum_function_uz; - int parser_instance_number_ux=0; - int parser_instance_number_uy=0; - int parser_instance_number_uz=0; + WarpXParser parser_ux; + WarpXParser parser_uy; + WarpXParser parser_uz; }; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 3c120a7c1..01d122ac6 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -74,17 +74,13 @@ ParseDensityProfile::ParseDensityProfile(std::string parse_density_function) { my_constants.ReadParameters(); parse_density_function = my_constants.replaceStringValue(parse_density_function); - const std::string s_var = "x,y,z"; - parser_instance_number = parser_initialize_function(parse_density_function.c_str(), - parse_density_function.length(), - s_var.c_str(), - s_var.length()); + parser_density.define(parse_density_function); + parser_density.registerVariables({"x","y","z"}); } Real ParseDensityProfile::getDensity(Real x, Real y, Real z) const { - std::array<amrex::Real, 3> list_var = {x,y,z}; - return parser_evaluate_function(list_var.data(), 3, parser_instance_number); + return parser_density.eval(x,y,z); } ConstantMomentumDistribution::ConstantMomentumDistribution(Real ux, @@ -141,31 +137,23 @@ ParseMomentumFunction::ParseMomentumFunction(std::string parse_momentum_function _parse_momentum_function_uy(parse_momentum_function_uy), _parse_momentum_function_uz(parse_momentum_function_uz) { - const std::string s_var = "x,y,z"; my_constants.ReadParameters(); parse_momentum_function_ux = my_constants.replaceStringValue(parse_momentum_function_ux); parse_momentum_function_uy = my_constants.replaceStringValue(parse_momentum_function_uy); parse_momentum_function_uz = my_constants.replaceStringValue(parse_momentum_function_uz); - parser_instance_number_ux = parser_initialize_function(parse_momentum_function_ux.c_str(), - parse_momentum_function_ux.length(), - s_var.c_str(), - s_var.length()); - parser_instance_number_uy = parser_initialize_function(parse_momentum_function_uy.c_str(), - parse_momentum_function_uy.length(), - s_var.c_str(), - s_var.length()); - parser_instance_number_uz = parser_initialize_function(parse_momentum_function_uz.c_str(), - parse_momentum_function_uz.length(), - s_var.c_str(), - s_var.length()); + parser_ux.define(parse_momentum_function_ux); + parser_uy.define(parse_momentum_function_uy); + parser_uz.define(parse_momentum_function_uz); + parser_ux.registerVariables({"x","y","z"}); + parser_uy.registerVariables({"x","y","z"}); + parser_uz.registerVariables({"x","y","z"}); } void ParseMomentumFunction::getMomentum(vec3& u, Real x, Real y, Real z) { - std::array<amrex::Real, 3> list_var = {x,y,z}; - u[0] = parser_evaluate_function(list_var.data(), 3, parser_instance_number_ux); - u[1] = parser_evaluate_function(list_var.data(), 3, parser_instance_number_uy); - u[2] = parser_evaluate_function(list_var.data(), 3, parser_instance_number_uz); + u[0] = parser_ux.eval(x,y,z); + u[1] = parser_uy.eval(x,y,z); + u[2] = parser_uz.eval(x,y,z); } RandomPosition::RandomPosition(int num_particles_per_cell): diff --git a/Source/Initialization/WarpX_parser.F90 b/Source/Initialization/WarpX_parser.F90 deleted file mode 100644 index 4b95c5ac7..000000000 --- a/Source/Initialization/WarpX_parser.F90 +++ /dev/null @@ -1,1365 +0,0 @@ -! MODULE mod_interpret -! Parser to read and evaluate math expressions - -module mod_interpret -use iso_c_binding -use amrex_fort_module, only : amrex_real -use amrex_error_module, only : amrex_error -implicit none -INTEGER, parameter :: plus =1, & - minus =2, & - multiply =3, & - divide =4, & - power =5, & - greaterthan=6, & - lessthan =7, & - exponential=8, & - logarithm =9, & - sine =10, & - cosine =11, & - tangent =12, & - square_root=13, & - arccos =14, & - arcsin =15, & - arctan =16, & - sinhyp =17, & - coshyp =18, & - tanhyp =19, & - logten =20 - -CHARACTER(5), DIMENSION(0:20) :: coper -INTEGER, parameter :: w_parenthesis = 1, & - w_number = 2, & - w_operator = 3, & - w_intrinsic = 4, & - w_variable = 5 -TYPE :: operation_type - INTEGER :: a,b,op -END TYPE operation_type - -TYPE :: res_type - INTEGER :: nb_op - REAL(amrex_real) :: value - TYPE(res_type), DIMENSION(:), POINTER :: res - LOGICAL, DIMENSION(:), POINTER :: l_res, l_res_do - TYPE(operation_type), DIMENSION(:), POINTER :: operation - LOGICAL :: a_res,a_l_res,a_l_res_do,a_operation -END TYPE res_type - -TYPE :: res_type_array - INTEGER :: nb_op - REAL(amrex_real), DIMENSION(:), POINTER :: value - TYPE(res_type_array), DIMENSION(:), POINTER :: res - LOGICAL, DIMENSION(:), POINTER :: l_res, l_res_do - TYPE(operation_type), DIMENSION(:), POINTER :: operation - LOGICAL :: a_res,a_l_res,a_l_res_do,a_operation -END TYPE res_type_array - -TYPE(res_type_array) :: res_array - -INTEGER :: nb_res=0 -TYPE(res_type), DIMENSION(10) :: table_of_res - -contains - -RECURSIVE subroutine del_res(resin) -implicit none -TYPE(res_type), INTENT(INOUT) :: resin -INTEGER :: j - IF(resin%a_res) then - do j = LBOUND(resin%res,1),UBOUND(resin%res,1) - call del_res(resin%res(j)) - end do - DEALLOCATE(resin%res) - resin%a_res=.FALSE. - END if - IF(resin%a_l_res) then - DEALLOCATE(resin%l_res) - resin%a_l_res=.FALSE. - END if - IF(resin%a_l_res_do) then - DEALLOCATE(resin%l_res_do) - resin%a_l_res_do=.FALSE. - END if - IF(resin%a_operation) then - DEALLOCATE(resin%operation) - resin%a_operation=.FALSE. - END if -return -end subroutine del_res - -RECURSIVE subroutine del_res_array(resin) -implicit none -TYPE(res_type_array), INTENT(INOUT) :: resin -INTEGER :: j - IF(resin%a_res) then - do j = LBOUND(resin%res,1),UBOUND(resin%res,1) - call del_res_array(resin%res(j)) - end do - DEALLOCATE(resin%res) - resin%a_res=.FALSE. - END if - IF(resin%a_l_res) then - DEALLOCATE(resin%l_res) - resin%a_l_res=.FALSE. - END if - IF(resin%a_l_res_do) then - DEALLOCATE(resin%l_res_do) - resin%a_l_res_do=.FALSE. - END if - IF(resin%a_operation) then - DEALLOCATE(resin%operation) - resin%a_operation=.FALSE. - END if - if (associated(resin%value)) then - nullify(resin%value) - end if -return -end subroutine del_res_array - -subroutine init_interpret() -implicit none - -coper(0 )='###' -coper(plus )='+' -coper(minus )='-' -coper(multiply )='*' -coper(divide )='/' -coper(power )='**' -coper(square_root)='sqrt' -coper(exponential)='exp' -coper(logarithm )='log' -coper(sine )='sin' -coper(cosine )='cos' -coper(tangent )='tan' -coper(arccos )='acos' -coper(arcsin )='asin' -coper(arctan )='atan' -coper(sinhyp )='sinh' -coper(coshyp )='cosh' -coper(tanhyp )='tanh' -coper(logten )='log10' -coper(greaterthan)='>' -coper(lessthan )='<' - -return -end subroutine init_interpret - -recursive function calc_res(res,list_var,root) RESULT(calc_res_r) -implicit none -REAL(amrex_real) :: calc_res_r -TYPE(res_type), INTENT(IN OUT) :: res -REAL(amrex_real), DIMENSION(:), OPTIONAL :: list_var -LOGICAL, INTENT(IN), OPTIONAL :: root - -INTEGER :: j,ia,ib,op -REAL(amrex_real) :: a,b - -do j = 1, res%nb_op - ia = res%operation(j)%a - ib = res%operation(j)%b - op = res%operation(j)%op - IF(op<0) then - IF(op>-100) res%res(ia)%value=list_var(-op) - else - IF(res%l_res(ia).and.res%l_res_do(ia)) then - IF(PRESENT(list_var)) then - a = calc_res(res=res%res(ia),list_var=list_var,root=.false.) - else - a = calc_res(res=res%res(ia),root=.false.) - END if - res%l_res_do(ia)=.false. - else - a = res%res(ia)%value - END if - IF(res%l_res(ib).and.res%l_res_do(ib)) then - IF(PRESENT(list_var)) then - b = calc_res(res=res%res(ib),list_var=list_var,root=.false.) - else - b = calc_res(res=res%res(ib),root=.false.) - END if - ! This might prevent OMP threading - res%l_res_do(ib)=.false. - else - b = res%res(ib)%value - END if - ! This might prevent OMP threading - res%res(ia)%value=eval(a,b,op) - END if -end do -IF(res%nb_op==0) then - ia=1 - calc_res_r = res%res(ia)%value -else - calc_res_r = res%res(ia)%value -END if - -IF(.not.PRESENT(root)) call calc_res_init(res) - -end function calc_res - -recursive subroutine calc_res_init(res) -implicit none -TYPE(res_type), INTENT(IN OUT) :: res - -INTEGER :: j,ia,ib,op - -do j = 1, res%nb_op - ia = res%operation(j)%a - ib = res%operation(j)%b - op = res%operation(j)%op - IF(op>=0) then - IF(res%l_res(ia).and..not.res%l_res_do(ia)) then - call calc_res_init(res%res(ia)) - res%l_res_do(ia)=.true. - END if - IF(res%l_res(ib).and..not.res%l_res_do(ib)) then - call calc_res_init(res%res(ib)) - res%l_res_do(ib)=.true. - END if - END if -end do - -end subroutine calc_res_init - -recursive function calc_res_array(res,list_values,root) RESULT(calc_res_array_r) -implicit none -REAL(amrex_real), dimension(:), pointer :: calc_res_array_r -TYPE(res_type_array), INTENT(IN OUT) :: res -REAL(amrex_real), DIMENSION(:,:), OPTIONAL :: list_values -LOGICAL, INTENT(IN), OPTIONAL :: root - -INTEGER :: j,ia,ib,op,n -REAL(amrex_real),dimension(:),pointer :: a,b - -do j = 1, res%nb_op - ia = res%operation(j)%a - ib = res%operation(j)%b - op = res%operation(j)%op - IF(op<0) then - IF(op>-100) then - n = size(list_values,2) - allocate(res%res(ia)%value(n)) - res%res(ia)%value=list_values(-op,:) - end if - else - IF(res%l_res(ia).and.res%l_res_do(ia)) then - IF(PRESENT(list_values)) then - a => calc_res_array(res=res%res(ia),list_values=list_values,root=.false.) - else - a => calc_res_array(res=res%res(ia),root=.false.) - END if - res%l_res_do(ia)=.false. - else - a => res%res(ia)%value - END if - IF(res%l_res(ib).and.res%l_res_do(ib)) then - IF(PRESENT(list_values)) then - b => calc_res_array(res=res%res(ib),list_values=list_values,root=.false.) - else - b => calc_res_array(res=res%res(ib),root=.false.) - END if - res%l_res_do(ib)=.false. - else - b => res%res(ib)%value - END if - res%res(ia)%value=>eval_array(a,b,op) - END if -end do - -IF(res%nb_op==0) then - ia=1 - calc_res_array_r => res%res(ia)%value -else - calc_res_array_r => res%res(ia)%value -END if - -IF(.not.PRESENT(root)) call calc_res_array_init(res) - -end function calc_res_array - -recursive subroutine calc_res_array_init(res) -implicit none -TYPE(res_type_array), INTENT(IN OUT) :: res - -INTEGER :: j,ia,ib,op - -do j = 1, res%nb_op - ia = res%operation(j)%a - ib = res%operation(j)%b - op = res%operation(j)%op - IF(op>=0) then - IF(res%l_res(ia).and..not.res%l_res_do(ia)) then - call calc_res_array_init(res%res(ia)) - res%l_res_do(ia)=.true. - END if - IF(res%l_res(ib).and..not.res%l_res_do(ib)) then - call calc_res_array_init(res%res(ib)) - res%l_res_do(ib)=.true. - END if - END if -end do - -end subroutine calc_res_array_init - -recursive subroutine eval_res(res,exprin,int_op,list_var,l_root) -implicit none -TYPE(res_type) :: res -CHARACTER(LEN=*), INTENT(IN) :: exprin -INTEGER, OPTIONAL, INTENT(IN) :: int_op -CHARACTER(*), DIMENSION(:),INTENT(IN), OPTIONAL :: list_var -LOGICAL, OPTIONAL :: l_root - -CHARACTER(LEN=LEN(exprin)) :: expr_old -CHARACTER(LEN=LEN(exprin)+2) :: expr -INTEGER, parameter :: nop_tot=25 -INTEGER :: ln, i, j,ja -INTEGER :: jtot,jres,jop -INTEGER, DIMENSION(0:nop_tot) :: next_l,next_r,op,what -LOGICAL, ALLOCATABLE, DIMENSION(:) :: l_res - -IF(.not.PRESENT(l_root)) call del_res(res) - -res%a_res=.false. -res%a_l_res=.false. -res%a_l_res_do=.false. -res%a_operation=.false. - -op=0 -what=0 - -expr_old = TRIM(ADJUSTL(C_replace(exprin,' ',''))) -IF(expr_old(1:1)=='-') then - expr='0'//expr_old -else - expr='0+'//expr_old -END if -ln = LEN(TRIM(ADJUSTL(expr))) -next_r(0) = 0 -i = 1 -j = 1 -jres = 0 -do WHILE(i<ln+1) - IF(PRESENT(list_var)) then - call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j),list_var) - else - call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j)) - END if - next_r(j) = next_r(j)+next_r(j-1) - next_l(j) = next_r(j-1)+1+next_l(j) - i = next_r(j)+1 - IF(what(j)/=w_operator) jres=jres+1 - j = j+1 -end do -jtot = j-1 -IF(PRESENT(int_op)) then - ALLOCATE(res%res(0:jres),res%operation(jres),res%l_res(0:jres),res%l_res_do(0:jres)) - res%a_res=.true. - res%a_l_res=.true. - res%a_l_res_do=.true. - res%a_operation=.true. - res%res(0:jres)%a_res=.false. - res%res(0:jres)%a_l_res=.false. - res%res(0:jres)%a_l_res_do=.false. - res%res(0:jres)%a_operation=.false. - res%l_res(0)=.false. - res%l_res_do=.true. - res%res(0)%value=0. - res%nb_op=jres - res%operation(jres)%op=int_op - res%operation(jres)%a=1 - res%operation(jres)%b=0 - res%l_res(jres)=.false. - nullify(res%res(0)%res) - nullify(res%res(0)%l_res) - nullify(res%res(0)%l_res_do) - nullify(res%res(0)%operation) -else - ALLOCATE(res%res(jres),res%operation(jres-1),res%l_res(jres),res%l_res_do(jres)) - res%a_res=.true. - res%a_l_res=.true. - res%a_l_res_do=.true. - res%a_operation=.true. - res%res(1:jres)%a_res=.false. - res%res(1:jres)%a_l_res=.false. - res%res(1:jres)%a_l_res_do=.false. - res%res(1:jres)%a_operation=.false. - res%nb_op=jres-1 - res%l_res_do=.true. -END if -ALLOCATE(l_res(jres)) -l_res(:) = .true. - -res%res(:)%value=0. - -jres = 0 -do j=1, jtot - select case (what(j)) - case (w_number) - jres = jres+1 - ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1)) - res%res(jres)%a_res=.true. - res%res(jres)%a_l_res=.true. - res%res(jres)%a_l_res_do=.true. - res%res(jres)%a_operation=.true. - res%res(jres)%res(1)%a_res=.false. - res%res(jres)%res(1)%a_l_res=.false. - res%res(jres)%res(1)%a_l_res_do=.false. - res%res(jres)%res(1)%a_operation=.false. - res%res(jres)%nb_op=1 - res%res(jres)%l_res(1) = .false. - res%res(jres)%l_res_do(1) = .true. - res%res(jres)%operation(1)%op = -100-op(j) - res%res(jres)%operation(1)%a = 1 - res%res(jres)%operation(1)%b = 0 - res%l_res(jres)=.true. - res%res(jres)%res(1)%value = char2real(expr(next_l(j):next_r(j))) - nullify(res%res(jres)%res(1)%res) - nullify(res%res(jres)%res(1)%l_res) - nullify(res%res(jres)%res(1)%l_res_do) - nullify(res%res(jres)%res(1)%operation) - case (w_parenthesis) - jres = jres+1 - IF(PRESENT(list_var)) then - call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),list_var=list_var,l_root=.false.) - else - call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),l_root=.false.) - endif - res%l_res(jres)=.true. - case (w_intrinsic) - jres = jres+1 - IF(PRESENT(list_var)) then - call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),list_var,l_root=.false.) - else - call eval_res(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),l_root=.false.) - END if - res%l_res(jres)=.true. - case (w_variable) - jres = jres+1 - ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1)) - res%res(jres)%a_res=.true. - res%res(jres)%a_l_res=.true. - res%res(jres)%a_l_res_do=.true. - res%res(jres)%a_operation=.true. - res%res(jres)%res(1)%a_res=.false. - res%res(jres)%res(1)%a_l_res=.false. - res%res(jres)%res(1)%a_l_res_do=.false. - res%res(jres)%res(1)%a_operation=.false. - res%res(jres)%nb_op=1 - res%res(jres)%l_res(1) = .false. - res%res(jres)%l_res_do(1) = .true. - res%res(jres)%operation(1)%op = -op(j) - res%res(jres)%operation(1)%a = 1 - res%res(jres)%operation(1)%b = 0 - res%l_res(jres)=.true. - nullify(res%res(jres)%res(1)%res) - nullify(res%res(jres)%res(1)%l_res) - nullify(res%res(jres)%res(1)%l_res_do) - nullify(res%res(jres)%res(1)%operation) - case (w_operator) - case default - end select -end do -jop = 1 -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==power) then - res%operation(jop)%op = power - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END IF - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==divide) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==multiply) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==plus.or.op(j)==minus) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==greaterthan.or.op(j)==lessthan) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do - -DEALLOCATE(l_res) - -return -end subroutine eval_res - -recursive subroutine eval_res_array(res,exprin,int_op,list_var,l_root) -implicit none -TYPE(res_type_array) :: res -CHARACTER(LEN=*), INTENT(IN) :: exprin -INTEGER, OPTIONAL, INTENT(IN) :: int_op -CHARACTER(*), DIMENSION(:),INTENT(IN), OPTIONAL :: list_var -LOGICAL, OPTIONAL :: l_root - -CHARACTER(LEN=LEN(exprin)) :: expr_old -CHARACTER(LEN=LEN(exprin)+2) :: expr -INTEGER, parameter :: nop_tot=25 -INTEGER :: ln, i, j,ja -INTEGER :: jtot,jres,jop -INTEGER, DIMENSION(0:nop_tot) :: next_l,next_r,op,what -LOGICAL, ALLOCATABLE, DIMENSION(:) :: l_res - -IF(.not.PRESENT(l_root)) call del_res_array(res) - -res%a_res=.false. -res%a_l_res=.false. -res%a_l_res_do=.false. -res%a_operation=.false. - -op=0 -what=0 - -expr_old = TRIM(ADJUSTL(C_replace(exprin,' ',''))) -IF(expr_old(1:1)=='-') then - expr='0'//expr_old -else - expr='0+'//expr_old -END if -ln = LEN(TRIM(ADJUSTL(expr))) -next_r(0) = 0 -i = 1 -j = 1 -jres = 0 -do WHILE(i<ln+1) - IF(PRESENT(list_var)) then - call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j),list_var) - else - call what_is_next(expr(i:ln),what(j),next_l(j),next_r(j),op(j)) - END if - next_r(j) = next_r(j)+next_r(j-1) - next_l(j) = next_r(j-1)+1+next_l(j) - i = next_r(j)+1 - IF(what(j)/=w_operator) jres=jres+1 - j = j+1 -end do -jtot = j-1 -IF(PRESENT(int_op)) then - ALLOCATE(res%res(0:jres),res%operation(jres),res%l_res(0:jres),res%l_res_do(0:jres)) - res%a_res=.true. - res%a_l_res=.true. - res%a_l_res_do=.true. - res%a_operation=.true. - res%res(0:jres)%a_res=.false. - res%res(0:jres)%a_l_res=.false. - res%res(0:jres)%a_l_res_do=.false. - res%res(0:jres)%a_operation=.false. - res%l_res(0)=.false. - res%l_res_do=.true. - nullify(res%res(0)%value) - res%nb_op=jres - res%operation(jres)%op=int_op - res%operation(jres)%a=1 - res%operation(jres)%b=0 - res%l_res(jres)=.false. - nullify(res%res(0)%res) - nullify(res%res(0)%l_res) - nullify(res%res(0)%l_res_do) - nullify(res%res(0)%operation) -else - ALLOCATE(res%res(jres),res%operation(jres-1),res%l_res(jres),res%l_res_do(jres)) - res%a_res=.true. - res%a_l_res=.true. - res%a_l_res_do=.true. - res%a_operation=.true. - res%res(1:jres)%a_res=.false. - res%res(1:jres)%a_l_res=.false. - res%res(1:jres)%a_l_res_do=.false. - res%res(1:jres)%a_operation=.false. - res%nb_op=jres-1 - res%l_res_do=.true. -END if -ALLOCATE(l_res(jres)) -l_res(:) = .true. - -jres = 0 -do j=1, jtot - select case (what(j)) - case (w_number) - jres = jres+1 - ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1)) - res%res(jres)%a_res=.true. - res%res(jres)%a_l_res=.true. - res%res(jres)%a_l_res_do=.true. - res%res(jres)%a_operation=.true. - res%res(jres)%res(1)%a_res=.false. - res%res(jres)%res(1)%a_l_res=.false. - res%res(jres)%res(1)%a_l_res_do=.false. - res%res(jres)%res(1)%a_operation=.false. - res%res(jres)%nb_op=1 - res%res(jres)%l_res(1) = .false. - res%res(jres)%l_res_do(1) = .true. - res%res(jres)%operation(1)%op = -100-op(j) - res%res(jres)%operation(1)%a = 1 - res%res(jres)%operation(1)%b = 0 - res%l_res(jres)=.true. - allocate(res%res(jres)%res(1)%value(1)) - res%res(jres)%res(1)%value = char2real(expr(next_l(j):next_r(j))) - nullify(res%res(jres)%res(1)%res) - nullify(res%res(jres)%res(1)%l_res) - nullify(res%res(jres)%res(1)%l_res_do) - nullify(res%res(jres)%res(1)%operation) - case (w_parenthesis) - jres = jres+1 - IF(PRESENT(list_var)) then - call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),list_var=list_var,l_root=.false.) - else - call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),l_root=.false.) - endif - res%l_res(jres)=.true. - case (w_intrinsic) - jres = jres+1 - IF(PRESENT(list_var)) then - call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),list_var,l_root=.false.) - else - call eval_res_array(res%res(jres),expr(next_l(j)+1:next_r(j)-1),op(j),l_root=.false.) - END if - res%l_res(jres)=.true. - case (w_variable) - jres = jres+1 - ALLOCATE(res%res(jres)%res(1),res%res(jres)%l_res(1),res%res(jres)%l_res_do(1),res%res(jres)%operation(1)) - res%res(jres)%a_res=.true. - res%res(jres)%a_l_res=.true. - res%res(jres)%a_l_res_do=.true. - res%res(jres)%a_operation=.true. - res%res(jres)%res(1)%a_res=.false. - res%res(jres)%res(1)%a_l_res=.false. - res%res(jres)%res(1)%a_l_res_do=.false. - res%res(jres)%res(1)%a_operation=.false. - res%res(jres)%nb_op=1 - res%res(jres)%l_res(1) = .false. - res%res(jres)%l_res_do(1) = .true. - res%res(jres)%operation(1)%op = -op(j) - res%res(jres)%operation(1)%a = 1 - res%res(jres)%operation(1)%b = 0 - res%l_res(jres)=.true. - nullify(res%res(jres)%res(1)%res) - nullify(res%res(jres)%res(1)%l_res) - nullify(res%res(jres)%res(1)%l_res_do) - nullify(res%res(jres)%res(1)%operation) - case (w_operator) - case default - end select -end do -jop = 1 -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==power) then - res%operation(jop)%op = power - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END IF - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==divide) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==multiply) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==plus.or.op(j)==minus) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do -do j=jtot,1,-1 - IF(what(j)==w_operator) then - IF(op(j)==greaterthan.or.op(j)==lessthan) then - res%operation(jop)%op = op(j) - ja=j/2 - do WHILE(.not.l_res(ja)) - ja=ja-1 - end do - res%operation(jop)%a = ja - res%operation(jop)%b = j/2+1 - l_res(j/2+1)=.false. - jop = jop+1 - END if - END if -end do - -DEALLOCATE(l_res) - -return -end subroutine eval_res_array - -subroutine what_is_next(exprin,what,next_l,next_r,op,list_var) -implicit none -CHARACTER(LEN=*), INTENT(IN) :: exprin -INTEGER, INTENT(OUT) :: what,next_l,next_r,op -CHARACTER(*), DIMENSION(:),INTENT(IN), optional :: list_var - -CHARACTER(LEN=LEN(exprin)) :: expr - -INTEGER :: ln, i, asc, istart,ln2,ierror,jv - -ln = LEN_TRIM(ADJUSTL(exprin)) -expr = TRIM(ADJUSTL(exprin)) -istart = LEN(exprin)-LEN(ADJUSTL(exprin)) -next_l = istart - -i = 1 -asc = IACHAR(expr(1:1)) -select case (asc) - ! parenthesis ( - case (40) - what = w_parenthesis - next_r = find_close_bracket(expr(i:ln)) - ! alphabet - case (65:90,97:122) - ln2 = ln-i+1 - op = find_intrinsic(expr(i:ln),ln2,ierror) - IF(ierror==0) then - what = w_intrinsic - IF( expr(i+ln2:i+ln2) .ne.'(' ) then - WRITE(*,*) 'Error in intrinsic ',exprin,': missing parenthesis.' - stop - END if - next_l = next_l+ln2 - next_r = ln2+find_close_bracket(expr(i+ln2:ln)) - ELSEIF(PRESENT(list_var)) then - what = w_variable - jv = 0 - op = 0 - do WHILE(jv+1<=SIZE(list_var)) - jv=jv+1 - IF(expr(i:i+ln2-1)==TRIM(list_var(jv)(:))) then - op=jv - END if - end do - IF(op==0) then - WRITE(*,*) 'Error: variable ',exprin(i:i+ln2-1),' not allowed.' - stop - END if - next_l = next_l+ln2 - next_r = ln2 - else - WRITE(*,*) 'Error: intrinsic ',exprin(i:i+ln2-1),' does not exists.' - stop - END if - ! number - case (48:57) - what = w_number - next_r = next_l+find_number_end(expr(i:ln)) - ! operator - case (42,43,45,47,60,62) - what = w_operator - op = find_operator(expr(i:i+1)) - IF(op==power) THEN - next_r = next_l+2 - else - next_r = next_l+1 - END if - case default -end select - -return -end subroutine what_is_next - -function find_close_bracket(a) -implicit none -INTEGER :: find_close_bracket -CHARACTER(*), INTENT(IN) :: a -INTEGER :: level, i, ln -level = 0 -ln=LEN(a) -do i = 1, ln - select case (a(i:i)) - case ('(') - level=level+1 - case (')') - level=level-1 - case default - end select - IF(level==0) exit -end do -IF(i>ln) then - WRITE(*,*) 'Error, no closing bracket in expression ',a - stop -END if - -find_close_bracket = i - -return -end function find_close_bracket - -function find_intrinsic(a,ln,ierror) -implicit none -INTEGER :: find_intrinsic -CHARACTER(*), INTENT(IN) :: a -CHARACTER(LEN=LEN(a)) :: ac -INTEGER, INTENT(IN OUT) :: ln -INTEGER, INTENT(OUT) :: ierror -INTEGER :: i - -ierror=0 -do i = 1, ln - select case (IACHAR(a(i:i))) - case (65:90,97:122) ! alphabet - case default - exit - end select -end do -IF(i>ln+1) then - WRITE(*,*) 'Error in expression ',a - stop -END if - -ln = i-1 -ac(1:ln) = a(1:ln) - -do i = 1, ln - IF(IACHAR(a(i:i))<97) ac(i:i)=ACHAR(IACHAR(a(i:i))+32) -END do - -select case (ac(1:ln)) - case ('sqrt') - find_intrinsic = square_root - case ('sin') - find_intrinsic = sine - case ('cos') - find_intrinsic = cosine - case ('tan') - find_intrinsic = tangent - case ('exp') - find_intrinsic = exponential - case ('log','ln') - find_intrinsic = logarithm - case ('log10') - find_intrinsic = logten - case ('asin') - find_intrinsic = arcsin - case ('acos') - find_intrinsic = arccos - case ('atan') - find_intrinsic = arctan - case ('sinh') - find_intrinsic = sinhyp - case ('cosh') - find_intrinsic = coshyp - case ('tanh') - find_intrinsic = tanhyp - case default - ierror=1 -end select -return -end function find_intrinsic - -function find_operator(a) -implicit none -INTEGER :: find_operator -CHARACTER(*), INTENT(IN) :: a - - select case (IACHAR(a(1:1))) - case (42) ! '*' - IF(IACHAR(a(2:2))==42) then - find_operator = power - else - find_operator = multiply - END if - case (43) ! '+' - find_operator = plus - case (45) ! '-' - find_operator = minus - case (47) ! '/' - find_operator = divide - case (60) ! '<' - find_operator = lessthan - case (62) ! '>' - find_operator = greaterthan - case default - WRITE(*,*) 'Error, operator ',a(1:1),' does not exist' - stop - end select - - return -end function find_operator - -function find_number_end(a) -implicit none -INTEGER :: find_number_end -CHARACTER(*), INTENT(IN) :: a -INTEGER :: i -LOGICAL :: l_point,l_minus,l_e - -l_point=.FALSE. -l_minus=.FALSE. -l_e=.FALSE. - -do i = 1, len(a) - select case (IACHAR(a(i:i))) - case (48:57) - l_minus = .true. - case (46) ! . - IF(l_point) exit - l_point = .true. - l_minus = .true. - case (45) ! - - IF(l_minus) exit - l_minus = .true. - case (69,101) ! e - IF(l_e) exit - l_e = .true. - l_minus = .false. - case default - exit - end select -end do - -find_number_end = i-1 - -return -end function find_number_end - -function char2real(a) -implicit none -REAL(amrex_real) :: char2real -CHARACTER(*), INTENT(IN) :: a - -READ(a(:),*) char2real - -return -end function char2real - -FUNCTION eval(a,b,op) -implicit none -REAL(amrex_real) :: eval -REAL(amrex_real), INTENT(IN) :: a, b -INTEGER, INTENT(IN) :: op - -select case (op) - case (plus) - eval=a+b - case (minus) - eval=a-b - case (multiply) - eval=a*b - case (divide) - eval=a/b - case (power) - eval=a**b - case (greaterthan) - eval=0. - if (a > b) eval=1. - case (lessthan) - eval=0. - if (a < b) eval=1. - case (square_root) - eval=SQRT(a) - case (exponential) - eval=EXP(a) - case (logarithm) - eval=LOG(a) - case (sine) - eval=SIN(a) - case (cosine) - eval=COS(a) - case (tangent) - eval=TAN(a) - case (logten) - eval=log10(a) - case (arcsin) - eval=asin(a) - case (arccos) - eval=acos(a) - case (arctan) - eval=atan(a) - case (sinhyp) - eval=sinh(a) - case (coshyp) - eval=cosh(a) - case (tanhyp) - eval=tanh(a) - case default - WRITE(*,*) 'Error in eval, wrong operator ',op - stop -end select - -return -END function eval - - -FUNCTION eval_array(ai,bi,op) -implicit none -REAL(amrex_real), dimension(:), pointer :: eval_array -REAL(amrex_real), INTENT(IN), dimension(:), pointer :: ai, bi -REAL(amrex_real), dimension(:), pointer :: a, b -INTEGER, INTENT(IN) :: op -integer::i1,i2,n - -i1 = size(ai) -i2 = size(bi) -n = maxval((/i1,i2/)) - -select case (op) - case (plus,minus,multiply,divide,power, greaterthan, lessthan) - - if (i1<n) then - allocate(a(n)) - a(:) = ai(1) - else - a => ai - end if - if (i2<n) then - allocate(b(n)) - b(:) = bi(1) - else - b => bi - end if - - case default - a => ai - -end select - -allocate(eval_array(n)) -eval_array=0. - -select case (op) - case (plus) - eval_array=a+b - case (minus) - eval_array=a-b - case (multiply) - eval_array=a*b - case (divide) - eval_array=a/b - case (power) - eval_array=a**b - case (square_root) - eval_array=SQRT(a) - case (exponential) - eval_array=EXP(a) - case (logarithm) - eval_array=LOG(a) - case (sine) - eval_array=SIN(a) - case (cosine) - eval_array=COS(a) - case (tangent) - eval_array=TAN(a) - case (logten) - eval_array=log10(a) - case (arcsin) - eval_array=asin(a) - case (arccos) - eval_array=acos(a) - case (arctan) - eval_array=atan(a) - case (sinhyp) - eval_array=sinh(a) - case (coshyp) - eval_array=cosh(a) - case (tanhyp) - eval_array=tanh(a) - case default - WRITE(*,*) 'Error in eval_array, wrong operator ',op - stop -end select - -return -END function eval_array - -function C_up2low(cline) -CHARACTER(*) :: cline -CHARACTER(LEN=LEN(cline)) :: C_up2low -INTEGER :: ln, i, j - ln = LEN(cline) - do i = 1, ln - j = IACHAR(cline(i:i)) - IF(j<91.and.j>64) then - C_up2low(i:i)=ACHAR(j+32) - else - C_up2low(i:i) = cline(i:i) - END if - END do -end function C_up2low - -RECURSIVE FUNCTION C_REPLACE(ST,R1,R2) RESULT(C_RES) -IMPLICIT NONE -CHARACTER(LEN=*), INTENT(IN) :: ST -CHARACTER(LEN=LEN(st)) :: c_res -CHARACTER(LEN=*), INTENT(IN) :: R1 -CHARACTER(LEN=*), INTENT(IN) :: R2 - -INTEGER :: i,ll - - c_res = ST - i=INDEX(TRIM(st),r1) - IF(i==0) return - ll = LEN(r2) - IF(ll>0) c_res(i:i+ll-1) = r2 - c_res(i+ll:) = c_replace(c_res(i+LEN(r1):),r1,r2) - - -RETURN -END FUNCTION C_REPLACE - -FUNCTION C_nboccur(ST,R1) -IMPLICIT NONE -CHARACTER(LEN=*), INTENT(IN) :: ST -integer :: c_nboccur -CHARACTER(LEN=*), INTENT(IN) :: R1 - -INTEGER :: i,is - - C_nboccur=0 - is=0 -10 i = INDEX(st(is+1:),r1) - IF(i/=0) then - C_nboccur = C_nboccur+1 - is=is+i+LEN(r1) - GO TO 10 - END if - -RETURN -END FUNCTION C_nboccur - -function stringlist(strin) RESULT(res) -CHARACTER(*), INTENT(IN) :: strin -INTEGER :: i1, nwords, i, i2 - -CHARACTER(LEN(strin)+1) :: str -CHARACTER(80), DIMENSION(:), pointer :: res - -str=ADJUSTL(strin)//' ' -i1=1 -nwords = 0 -do WHILE(TRIM(str(i1:))/='') - i2 = 1 - i2 = INDEX(str(i1:),' ') - do WHILE(i2==1) - i1 = i1 + 1 - i2 = INDEX(str(i1:),' ') - end do - i1=i1+i2 - nwords = nwords + 1 -end do -ALLOCATE(res(nwords)) -res(:) = '' -i1=1 -do i = 1, nwords - i2 = INDEX(str(i1:),' ') - do WHILE(i2==1) - i1 = i1 + 1 - i2 = INDEX(str(i1:),' ') - end do - res(i)(1:i2-1) = str(i1:i1+i2-1) - i1=i1+i2 -end do - -return -end function stringlist - -SUBROUTINE csv2list(strin, strinlen, lofstr) -CHARACTER(*), INTENT(IN) :: strin -INTEGER, INTENT(IN) :: strinlen -INTEGER :: isep, isep_rel, nwords, i, j -INTEGER, DIMENSION(strinlen) :: word_start, word_end -CHARACTER(strinlen) :: str -CHARACTER(len=10), DIMENSION(:), allocatable, INTENT(inout) :: lofstr -nwords = 1 -str = strin -isep_rel = INDEX(str, ',') -isep = isep_rel -word_start(nwords) = 1 -DO WHILE(isep_rel > 0) - word_end(nwords) = isep-1 - nwords = nwords + 1 - word_start(nwords) = isep+1 - str = strin(isep+1:strinlen) - isep_rel = INDEX(str, ',') - isep = isep_rel + isep -ENDDO -word_end(nwords) = strinlen - -ALLOCATE( lofstr(nwords) ) - -DO i=1, nwords - lofstr(i) = "" - DO j = 1, word_end(i)-word_start(i)+1 - lofstr(i)(j:j) = strin(word_start(i)+j-1:word_start(i)+j-1) - end DO -ENDDO -END SUBROUTINE csv2list -end module mod_interpret - -! --------------------- -! MODULE parser_wrapper -! Wrappers to use the parser defined in MODULE mod_interpret. -MODULE parser_wrapper -USE iso_c_binding -USE amrex_fort_module, only : amrex_real -USE amrex_error_module, only : amrex_error -USE mod_interpret, only : nb_res, table_of_res, calc_res, eval_res, csv2list - -IMPLICIT NONE - -CONTAINS - -! FUNCTION parser_initialize_function RESULT my_index_res -! Initialize a res_type and assign it a unique identifier. -! INPUTS: -!> instr_func : CHAR* for a mathematical expression -!> e.g. instr_func = "3*cos(x)+y". -!> instr_var : CHAR* for a comma-separated list of variables in instr_func -!> e.g. instr_var = "x,y". -! OUTPUT: -!> my_index_res : INT parser index. Necessary if this module is used for several -!> parsers (e.g. one for the electron density, one for the ion -!> density, one for the laser profile etc.). -FUNCTION parser_initialize_function(instr_func, length_func, instr_var, length_var) RESULT(my_index_res) & - bind(c,name='parser_initialize_function') - INTEGER(c_int), INTENT(IN), VALUE :: length_func, length_var - CHARACTER(kind=c_char), INTENT(IN) :: instr_func(length_func) - CHARACTER(kind=c_char), INTENT(IN) :: instr_var(length_var) - INTEGER :: i - CHARACTER(len=10), DIMENSION(:), allocatable :: lofstr - CHARACTER(len=length_func) :: str_func - CHARACTER(len=length_var) :: str_var - REAL(amrex_real), DIMENSION(1) :: list_var - INTEGER :: my_index_res - - nb_res = nb_res + 1 - my_index_res = nb_res - IF (nb_res>10) THEN - call amrex_error('Parser error: cannot have more than 10 parsers.') - ENDIF - - str_func = "" - DO i=1, length_func - str_func(i:i) = instr_func(i) - ENDDO - - str_var = "" - DO i=1, length_var - str_var(i:i) = instr_var(i) - ENDDO - - ! Convert variable list from csv string to list of strings. - CALL csv2list(str_var, len_trim(str_var), lofstr) - - ! Initialize the res object - CALL eval_res(table_of_res(my_index_res), str_func, list_var=lofstr) - -END FUNCTION parser_initialize_function - -! FUNCTION parser_evaluate_function RESULT out -! Evaluate parsed function -! INPUTS: -!> list_var : REAL* array of values for variables -!> e.g. list_var = (/3.14_8,2._8/). -!> nvar : INT number of variables -!> e.g. nvar = 2. -!> my_index_res : INT index of the res_type object in table table_of_res. -! OUTPUT: -!> out : REAL Result. -FUNCTION parser_evaluate_function(list_var, nvar, my_index_res) result(out) & - bind(c,name='parser_evaluate_function') - INTEGER, VALUE, INTENT(IN) :: nvar, my_index_res - REAL(amrex_real), INTENT(IN) :: list_var(1:nvar) - REAL(amrex_real) :: out - ! Evaluate parsed function in table_of_res(my_index_res), of type res_type - out = calc_res(table_of_res(my_index_res),list_var) - -END FUNCTION parser_evaluate_function - -END MODULE parser_wrapper |