module testfunctions_m ! Test functions use screen_io_m, only : get_stderr use precision_m, only : FP use constants_m, only : pi use equilibrium_m, only : equilibrium_t use circular_equilibrium_m, only : circular_t use numerical_equilibrium_m, only : numerical_t use slab_equilibrium_m, only : slab_t implicit none public :: testfun_u public :: testfun_dudrhon public :: testfun_co public :: testfun_lambda public :: testfun_xi public :: testfun_helm_u private :: numcomp_helm_u contains real(FP) function testfun_u(equi, x, y, phi) !! Test_function u class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: rho, rhon, theta, xn, yn select type(equi) type is(circular_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x, y) testfun_u = cos( 3.0_FP/2.0_FP*pi* rhon) * sin(4_FP*theta) + 1.3_FP type is (numerical_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) testfun_u = cos( 3.0_FP/2.0_FP*pi* rhon)+ 1.3_FP type is (slab_t) xn = (x-equi%xmin)/(equi%xmax-equi%xmin) yn = (y-equi%ymin)/(equi%ymax-equi%ymin) testfun_u = -((1-xn**2+Cos((5*Pi*xn)/2.D0))*Sin(1.2D0-2*Pi*yn)) class default write(get_stderr(), *) & 'error(testfun_u): equilibrium not valid as testcase' error stop end select end function real(FP) function testfun_dudrhon(equi, x, y, phi) !! Gradient normal to flux surface label, via finite difference on analytic solution with small grid distance class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP), parameter :: eps = 1.0E-8_FP real(FP) :: dudx, dudy, eradx, erady dudx = ( testfun_u(equi, x+eps, y , phi) - testfun_u(equi, x-eps, y , phi) ) / (2.0_FP*eps) dudy = ( testfun_u(equi, x, y+eps, phi) - testfun_u(equi, x, y-eps, phi) ) / (2.0_FP*eps) call equi%erad(x, y, phi, eradx, erady) testfun_dudrhon = dudx*eradx + dudy*erady end function real(FP) function testfun_co(equi, x, y, phi) !! Test function for polarisation coefficeint class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: rho, rhon, theta, xn, yn select type(equi) type is(circular_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x, y) testfun_co = 1.1_FP + cos(pi/2.0_FP * rhon) * cos(3_FP*theta) type is (numerical_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) testfun_co = 1.1_FP + cos(pi/2.0_FP * rhon) type is (slab_t) xn = (x-equi%xmin)/(equi%xmax-equi%xmin) yn = (y-equi%ymin)/(equi%ymax-equi%ymin) testfun_co = 1.1D0-(Cos(Pi*xn)*Sin(2.1D0-4*Pi*yn))/exp(xn) class default write(get_stderr(), *) & 'error(testcoeff_co): equilibrium not valid as testcase' error stop end select end function real(FP) function testfun_lambda(equi, x, y, phi) !! Test function for lambda class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: rho, rhon, theta, xn, yn select type(equi) type is(circular_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x, y) testfun_lambda = rho*sin(theta) type is (numerical_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) testfun_lambda = rho type is (slab_t) xn = (x-equi%xmin)/(equi%xmax-equi%xmin) yn = (y-equi%ymin)/(equi%ymax-equi%ymin) testfun_lambda = xn**2*(1-Sin(6.D-1-4*Pi*yn)) class default write(get_stderr(), *) & 'error(testfun_lambda): equilibrium not valid as testcase' error stop end select end function real(FP) function testfun_xi(equi, x, y, phi) !! Test function for xi class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: rho, rhon, theta, xn, yn select type(equi) type is(circular_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x,y) testfun_xi = sqrt(rho)*(2.0_FP+cos(theta)) type is (numerical_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) testfun_xi = sqrt(rho)*(2.0_FP) type is (slab_t) xn = (x-equi%xmin)/(equi%xmax-equi%xmin) yn = (y-equi%ymin)/(equi%ymax-equi%ymin) testfun_xi = 1+xn**3/2.D0 class default write(get_stderr(), *) & 'error(testfun_xi): equilibrium not valid as testcase' error stop end select end function real(FP) function testfun_helm_u(equi, x, y, phi) ! Analytic result of Helmholtz operator on u implicit none class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: rho, rhon, theta, xn, yn, nrmfc, r, t, rmin, rmax select type(equi) type is(circular_t) rho = equi%rho(x, y, phi) rhon = (rho-equi%rhomin) / (equi%rhomax - equi%rhomin) theta = equi%theta(x, y) r = rho t = theta rmin = equi%rhomin rmax = equi%rhomax testfun_helm_u = r*Sin(t)*(1.3D0+Cos((3*Pi*(r-rmin))/(2.D0*(rmax-rmin)))*Sin(4*t))-((2+Cos(t))*(-12*Cos((Pi*(r- & rmin))/(2.D0*(rmax-rmin)))*Cos((3*Pi*(r-rmin))/(2.D0*(rmax-rmin)))*Cos(4*t)*Sin(3*t)-16*Cos((3*Pi*(r & -rmin))/(2.D0*(rmax-rmin)))*(1.1D0+Cos((Pi*(r-rmin))/(2.D0*(rmax- & rmin)))*Cos(3*t))*Sin(4*t)))/r**1.5D0-((2+Cos(t))*((-9*Pi**2*r*Cos((3*Pi*(r-rmin))/(2.D0*(rmax- & rmin)))*(1.1D0+Cos((Pi*(r-rmin))/(2.D0*(rmax-rmin)))*Cos(3*t))*Sin(4*t))/(4.D0*(rmax-rmin)**2)- & (3*Pi*(1.1D0+Cos((Pi*(r-rmin))/(2.D0*(rmax-rmin)))*Cos(3*t))*Sin((3*Pi*(r-rmin))/(2.D0*(rmax- & rmin)))*Sin(4*t))/(2.D0*(rmax-rmin))+(3*Pi**2*r*Cos(3*t)*Sin((Pi*(r-rmin))/(2.D0*(rmax- & rmin)))*Sin((3*Pi*(r-rmin))/(2.D0*(rmax-rmin)))*Sin(4*t))/(4.D0*(rmax-rmin)**2)))/Sqrt(r) type is (numerical_t) testfun_helm_u = numcomp_helm_u(equi, x, y, phi) type is (slab_t) xn = (x-equi%xmin)/(equi%xmax-equi%xmin) yn = (y-equi%ymin)/(equi%ymax-equi%ymin) nrmfc = 1.0_FP/(equi%xmax-equi%xmin) testfun_helm_u = -(xn**2*(1-xn**2+Cos((5*Pi*xn)/2.D0))*(1-Sin(6.D-1-4*Pi*yn))*Sin(1.2D0-2*Pi*yn))-nrmfc**2*(1+ & xn**3/2.D0)*((8*Pi**2*Cos(Pi*xn)*(1-xn**2+Cos((5*Pi*xn)/2.D0))*Cos(2.1D0-4*Pi*yn)*Cos(1.2D0- & 2*Pi*yn))/exp(xn)+4*Pi**2*(1-xn**2+Cos((5*Pi*xn)/2.D0))*(1.1D0-(Cos(Pi*xn)*Sin(2.1D0- & 4*Pi*yn))/exp(xn))*Sin(1.2D0-2*Pi*yn)-(-2-(25*Pi**2*Cos((5*Pi*xn)/2.D0))/4.D0)*(1.1D0- & (Cos(Pi*xn)*Sin(2.1D0-4*Pi*yn))/exp(xn))*Sin(1.2D0-2*Pi*yn)-(-2*xn- & (5*Pi*Sin((5*Pi*xn)/2.D0))/2.D0)*((Cos(Pi*xn)*Sin(2.1D0-4*Pi*yn))/exp(xn)+(Pi*Sin(Pi*xn)*Sin(2.1D0- & 4*Pi*yn))/exp(xn))*Sin(1.2D0-2*Pi*yn)) class default write(get_stderr(), *) & 'error(testfun_helm_u): equilibrium not valid as testcase' error stop end select end function function numcomp_helm_u(equi, x, y, phi) result(res) !! Result of Helmholtz operator computed via finite difference at high accuracy class(equilibrium_t) :: equi real(FP), intent(in) :: x, y, phi real(FP) :: res real(FP), parameter :: eps = 1.0E-6_FP real(FP), parameter :: two_eps_sqrd = 2.0_FP*eps**2 real(FP) :: xlft, xrgt, ybtm, ytop real(FP) :: u_slf, u_lft, u_rgt, u_btm, u_top real(FP) :: co_slf, co_lft, co_rgt, co_btm, co_top real(FP) :: xi, lambda xlft = x - eps xrgt = x + eps ybtm = y - eps ytop = y + eps u_slf = testfun_u(equi, x, y , phi) u_lft = testfun_u(equi, xlft, y , phi) u_rgt = testfun_u(equi, xrgt, y , phi) u_btm = testfun_u(equi, x, ybtm, phi) u_top = testfun_u(equi, x, ytop, phi) co_slf = testfun_co(equi, x, y , phi) co_lft = testfun_co(equi, xlft, y , phi) + co_slf co_rgt = testfun_co(equi, xrgt, y , phi) + co_slf co_btm = testfun_co(equi, x, ybtm, phi) + co_slf co_top = testfun_co(equi, x, ytop, phi) + co_slf res = ( (co_btm + co_lft + co_rgt + co_top) * u_slf & - co_btm * u_btm & - co_lft * u_lft & - co_rgt * u_rgt & - co_top * u_top ) & / (two_eps_sqrd) xi = testfun_xi(equi, x, y, phi) lambda = testfun_lambda(equi, x, y, phi) res = xi*res + lambda*u_slf end function end module