Commit 17994b47 authored by Juan Caneses Marin (nfc)'s avatar Juan Caneses Marin (nfc)
Browse files

Used INTENT(IN/OUT/INOUT) in the definition and use of Interp1 and diff1, but... 

Used INTENT(IN/OUT/INOUT) in the definition and use of Interp1 and diff1, but no speedup yet with multiple threads
parent 893ec16d

InputFiles/xp_IonSlowingDown.in

+3 −3
Original line number Diff line number Diff line
@@ -12,12 +12,12 @@ in%nz = 501, ! N 


! Simulation conditions:

! =====================

in%Nparts = 50000,                                                              ! Total number of particles

in%Nparts = 500000,                                                              ! Total number of particles

in%Nsteps = 5000,                                                               ! Total number of time steps

in%dt     = 0.5E-7,                                                             ! Time step in [s]

in%zmax = +5.0,

in%zmin = -5.0,

in%threads_request = 1,

in%threads_request = 4,

in%padding = 4,



in%iSave = .true.,
@@ -31,7 +31,7 @@ in%iDrag = .false., ! 
! ==============

in%jstart = 1,                                                                  ! start frame

in%jend   = 5000,                                                               ! end frame

in%jincr  = 100,                                                                 ! increment

in%jincr  = 50,                                                                 ! increment



! Collision operator conditions:

! ==============================

src/Main.f90

+16 −2
Original line number Diff line number Diff line
@@ -127,9 +127,17 @@ WRITE(*,*) '' 


! Allocate memory to splines:

! ===========================================================================

! InitSpline takes in a variable of splType and performs the following:

! - Populates it scalar fields

! - Allocates memory for the profile data "y", indepedent coordinate "x", etc

CALL InitSpline(spline_Bz  ,in%nz,0._8,0._8,1,0._8)

CALL InitSpline(spline_ddBz,in%nz,0._8,0._8,1,0._8)

CALL InitSpline(spline_Phi ,in%nz,0._8,0._8,1,0._8)



! Allocate memory for spline_test variable:

! ==============================================================================

! Spline_test is a variable that holds up to 6 different test profiles "y" and 

! one "x" coordinate.

CALL InitSplineTest(spline_Test,in%nz)



! Allocate memory to main simulation variables:
@@ -140,9 +148,11 @@ ALLOCATE(ecount1(in%Nsteps),ecount2(in%Nsteps),ecount3(in%Nsteps),ecount4(in%Nst 


! Allocate memory to output variables:

! ==============================================================================

! Determine size of temporal snapshots to record:

jsize = (in%jend-in%jstart+1)/in%jincr

ALLOCATE(jrng(jsize))

ALLOCATE(zp_hist(in%Nparts,jsize),kep_hist(in%Nparts,jsize),xip_hist(in%Nparts,jsize),t_hist(jsize))



! Create array with the indices of the time steps to save:

jrng = (/ (j, j=in%jstart, in%jend, in%jincr) /)
@@ -157,7 +167,11 @@ ALLOCATE(fcurr(in%Nparts),fnew(in%Nparts)) 
fileName = trim(adjustl(in%BFieldFile))

fileName = trim(adjustl(in%BFieldFileDir))//fileName

fileName = trim(adjustl(in%rootDir))//fileName



! Populate profile data based on external file:

CALL ReadSpline(spline_Bz,fileName)



! Based on profile data. compute spline data:

CALL ComputeSpline(spline_Bz)



! Second derivative of the magnetic field:
@@ -165,7 +179,7 @@ spline_ddBz%x = spline_Bz%x 
spline_ddBz%y = spline_Bz%yp

CALL ComputeSpline(spline_ddBz)



! Predefined electric potential:

! Electric potential:

spline_Phi%x = spline_Bz%x

spline_Phi%y = 0

if (in%iPotential) then
@@ -356,7 +370,7 @@ TimeStepping: do j = 1,in%Nsteps 
    ! =====================================================================

    id = OMP_GET_THREAD_NUM()

    if (id .EQ. 0) then

      if (j .EQ. 1) then

      if (j .EQ. 50) then

	       oend_estimate = OMP_GET_WTIME()

         WRITE(*,*) 'Estimated compute time: ', in%Nsteps*(oend_estimate-ostart)/j,' [s]'

      end if

src/MoveParticlePack.f90

+54 −29
Original line number Diff line number Diff line
@@ -7,14 +7,17 @@ USE PhysicalConstants 
USE dataTYP



IMPLICIT NONE

! Define local variables

TYPE(splTYP) :: spline0, spline1

TYPE(inTYP)  :: in0

REAL(r8) :: zp0, kep0, xip0                  ! Position, kinetic energy and pitch of the ith particle



! Define type for interface arguments 

REAL(r8), INTENT(INOUT) :: zp0, kep0, xip0

TYPE(inTYP), INTENT(IN)  :: in0

TYPE(splTYP), INTENT(IN) :: spline0, spline1



! Local variables:

REAL(r8) :: zpnew, Xipnew, uparnew, upernew, munew  ! Position, kinetic energy and pitch of the ith particle

REAL(r8) :: m_t, dt



! Storage for the MoverParticle and RHS subroutines

! Storage for the MoverParticle and RHS subroutines:

REAL(r8) :: zp1, zp2, zp3

REAL(r8) :: K1, K2, K3, K4

REAL(r8) :: upar0, upar1, upar2, upar3
@@ -23,8 +26,8 @@ REAL(r8) :: mu0, mu1, mu2, mu3 
REAL(r8) :: M1, M2, M3, M4

REAL(r8) :: u2

REAL(r8) :: B, Phi

!REAL(r8) :: Interp1

REAL(r8) :: curv2

REAL(r8) :: Interp1

!REAL(r8) :: curv2



! Time step:

dt = in0%dt
@@ -39,8 +42,10 @@ upar0 = xip0*sqrt(2.*e_c*kep0/m_t) 
u2 = 2.*e_c*kep0/m_t



! Calculate initial magnetic moment:

!B = Interp1(zp0,spline0)

B = curv2(zp0,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

B = Interp1(zp0,spline0)

!B = curv2(zp0,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

!B = 1.



mu0 = 0.5*m_t*u2*(1 - xip0*xip0)/B



! Begin assembling RK4 solution:
@@ -65,8 +70,9 @@ uparnew = upar0 + ( (L1 + (2.*L2) + (2.*L3) + L4)/6. )*dt 
munew   = mu0 +   ( (M1 + (2.*M2) + (2.*M3) + M4)/6. )*dt



! Calculate the magnetic field at zpnew:

!B = Interp1(zpnew,spline0)

B = curv2(zpnew,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

B = Interp1(zpnew,spline0)

!B = curv2(zpnew,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

!B = 1.



! Based on new B and new mu, calculate new uper:

upernew = sqrt(2.*munew*B/m_t)
@@ -77,11 +83,12 @@ u2 = uparnew**2. + upernew**2. 
! New pitch angle due to new upar and new u:

Xipnew = uparnew/sqrt(u2)



! Output data:

zp0  = zpnew

xip0 = Xipnew

kep0 = 0.5*m_t*u2/e_c



return

RETURN

END SUBROUTINE MoveParticle



! =======================================================================================================
@@ -93,14 +100,18 @@ USE PhysicalConstants 
USE dataTYP



IMPLICIT NONE

! Define local variables:

TYPE(splTYP) :: spline0, spline1

TYPE(inTYP) :: in0

REAL(r8) :: zp0, upar0, mu0, K, L, M



! Define type of interface arguments:

REAL(r8), INTENT(IN) :: zp0, upar0, mu0

REAL(r8), INTENT(OUT) :: K, L, M

TYPE(inTYP), INTENT(IN) :: in0

TYPE(splTYP), INTENT(IN) :: spline0, spline1



! Local variables:

REAL(r8) :: dB, dPhi

REAL(r8) :: m_t, q_t

!REAL(r8) :: diff1

REAL(r8) :: curvd

REAL(r8) :: diff1

!REAL(r8) :: curvd



! Test particle mass:

m_t = in0%m_t
@@ -109,19 +120,21 @@ m_t = in0%m_t 
q_t = in0%q_t



! Calculate the magnetic field gradient at zp0:

!dB = diff1(zp0,spline0)

dB = curvd(zp0,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

dB = diff1(zp0,spline0)

!dB = curvd(zp0,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

!dB = .01



! Calculate electric potential gradient at zp0:

!dPhi = diff1(zp0,spline1)

dPhi = curvd(zp0,spline1%n,spline1%x,spline1%y,spline1%yp,spline1%sigma)

dPhi = diff1(zp0,spline1)

!dPhi = curvd(zp0,spline1%n,spline1%x,spline1%y,spline1%yp,spline1%sigma)

!dPhi = 0.



! Assign values to output variables:

K = upar0

L = -(1/m_t)*(mu0*dB + q_t*dPhi)

M = 0



return

RETURN

END SUBROUTINE RightHandSide



! =======================================================================================================
@@ -444,9 +457,15 @@ REAL(r8) FUNCTION Interp1(xq, spline0) 
  USE spline_fits



  IMPLICIT NONE

  TYPE(splTYP) :: spline0

  REAL(r8) :: xq, curv2

  

  ! Define type interface arguments:

  REAL(r8), INTENT(IN) :: xq

  TYPE(splTYP), INTENT(IN) :: spline0

  

  ! Local variables:

  REAL(r8) :: curv2



  ! Output data:

  Interp1 = curv2(xq,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

  !Interp1 = 1.
@@ -460,9 +479,15 @@ REAL(r8) FUNCTION diff1(xq, spline0) 
  USE spline_fits



  IMPLICIT NONE

  TYPE(splTYP) :: spline0

  REAL(r8) :: xq, curvd



  ! Define type interface arguments:

  REAL(r8), INTENT(IN) :: xq

  TYPE(splTYP), INTENT(IN) :: spline0



  ! Local variables:

  REAL(r8) :: curvd



  ! Output data:

  diff1 = curvd(xq,spline0%n,spline0%x,spline0%y,spline0%yp,spline0%sigma)

  ! diff1 = 0.1

src/compile.sh

deleted100644 → 0
+0 −1
Original line number Diff line number Diff line 
make -f Makefile_openMP.f