V3FIT
poincare.f
1  PROGRAM poincare
2 !
3 !=====================================================================================================================
4 !
5 ! Program = Poincare
6 ! Version = 1.10 03/30/10
7 ! 1.20 05/10/13 (add mpi)
8 !
9 ! Authors: Don A. Spong (original) and Raul Sanchez (updated), Fusion Energy Division, Oak Ridge National Laboratory
10 !
11 ! For problems, please contact: sanchezferlr@ornl.gov or spongda@ornl.gov
12 !
13 ! Description:
14 !
15 ! This program computes Poincare plots for a 3D magnetic field provided using flux-coordinates (s, theta, phi).
16 ! It is assumed that s in [0,1], u(=pol) in [0,2pi) and v(=tor) in [0,2pi). Values at u=2pi or v=2pi should
17 ! not be provided, since they are automatically set to the values at u=0 and v=0 to ensure periodicity
18 ! inside the program. s and u may be arbitrary flux coordinates, but it is assumed that v= geom. PHI, as in VMEC.
19 !
20 ! Poincare plots can be done on poloidal or toroidal cross-section, depending on input TYPE on command line,
21 ! which can be set to "poloidal" (theta is the independent variable) or "toroidal" (phi is the independent variable).
22 ! Choice should be done depending on whether Bu or Bv go through a zero. I.e., if BV is nonzero, choose "toroidal".
23 ! Angle for cross section can also be chosen on command line via the ANGLE input variable.
24 !
25 ! ** INPUT:The field information must be provided in an ASCII file whose first line contains the size of the 3D grid:
26 ! NV = No. of toroidal points; NS = No. of radial points; NU = No. of poloidal points
27 ! Then, data must follow in six colums, giving the cylindrical coordinates info:
28 ! R(s,u,v) Z(s,u,v) Phi(s,u,v) Br(s,u,v) Bz(s,u,v) BPhi(s,u,v)
29 ! The order assumed is that u is the fastest running, and v is the slowest running coordinate in the file.
30 !
31 ! ** EXECUTION: The program is run from the command line invoking:
32 ! xpoincare FILENAME.dat N.Toroidal transits N.surfaces N.orbits/surface N.steps/transit
33 ! Innermost Surface Outermost Surface TYPE ANGLE
34 !
35 ! Here, FILENAME.dat = name of file containing 3D field info (like the one produced by SIESTA code),
36 ! N. Toroidal transits = No. of 2*pi toroidal periods that must be followed per orbit
37 ! N. surfaces = No. of different s-values where orbits will be started
38 ! N. orbits/surface = No. of orbits to initialize, uniformly distributed, on each surface
39 ! N. steps/transit = the numerical integrator will integrate each transits in steps of size = 2*pi/(N.steps/transit)
40 ! Surfaces will be distributed between [Innermost Surface, Outermost Surface]
41 ! TYPE = "poloidal" (theta is ind.var.; choose if Bv nonzero) or "toroidal" (phi is ind. var.; choose if Bu nonzero)
42 ! ANGLE = angular cross-section where puncture plot computed. It is poloidal or toroidal depending on TYPE.
43 !
44 ! ** OUTPUT: The output is a 4-column ASCII file containing a list of the following values:
45 ! Theta s R Z
46 ! of all the crossings of the orbits with the plane PHI = ANGLE if TYPE = "toroidal".
47 ! Phi s R Phi
48 ! of all the crossings of the orbits with the plane THETA = ANGLE if TYPE = "poloidal".
49 ! The name of the file is simply FILENAME_puncture_TYPE.dat
50 !
51 ! ** HISTORY:
52 !
53 ! 10/24/09: First version released, based on a prior version by Don Spong. (Version: 1.00)
54 ! 03/30/10: Several minor bugs fixed. (Version: 1.10)
55 !
56 !=====================================================================================================================
57 !
58  USE stel_kinds
59  USE stel_constants
60  USE bspline ! Module in LIBSPLINE
61  USE spline3d_fit_coefs
62  USE lsode_quantities
63  USE mpi_params
64 #if defined(MPI_OPT)
65  USE mpi_inc
66 #endif
67  IMPLICIT NONE
68 
69  INTEGER:: i, j, k, iargc, numargs, isize, it,
70  & steps, num_orbs, icount, nsurfs, mn, scndry_steps
71  REAL(rprec):: s_test, theta_test, phi_test, delta,
72  & ang_in, ang_out, rdum, tstart, tend, s1, s2
73  REAL(rprec):: angle, phi_puncture, theta_puncture
74  REAL(rprec), PARAMETER :: relerr=1.d-8, abserr=1.d-8,
75  & eps=1.d-10
76  REAL(rprec),ALLOCATABLE, DIMENSION(:,:):: start_points
77  REAL(rprec),ALLOCATABLE, DIMENSION(:):: theta_nb_tot,
78  & phi_nb_tot, s_nb_tot, r_s_nb_tot, z_s_nb_tot
79  LOGICAL:: l_ex
80  INTEGER:: i_truncate, ilen, iargs
81  INTEGER :: local_num_orbs, norbs_tot, istat
82  CHARACTER(LEN=8):: type_in ! ="poloidal", THETA is ind. variable; = "toroidal", PHI is ind. variable
83  CHARACTER(LEN=100):: outfile
84  CHARACTER(LEN=60):: arg1, arg2, arg3, arg4, arg5,
85  & arg6, arg7, arg8, arg9, arg10
86  CHARACTER(LEN=13), PARAMETER:: banner = 'Version 1.20',
87  & update = 'May 10, 2013'
88  EXTERNAL rt_hand_side_dphi, rt_hand_side_dtheta, jacobian ! Used by the integrator LSODE
89 
90  CALL cpu_time(tstart)
91 
92 C
93 C Get NPES and myid. Requires initialization of MPI.
94 C
95 #if defined(MPI_OPT)
96  CALL mpi_init(ierr_mpi)
97  IF (ierr_mpi .NE. 0) THEN
98  WRITE(*,'("MPI_INIT returned IER =")') ierr_mpi
99  stop
100  ENDIF
101  CALL mpi_comm_size(mpi_comm_world, numprocs, ierr_mpi)
102  IF (ierr_mpi .NE. 0) THEN
103  WRITE(*,'("MPI_COMM_SIZE returned IER =")') ierr_mpi
104  stop
105  ENDIF
106  CALL mpi_comm_rank(mpi_comm_world, myid, ierr_mpi)
107  IF (ierr_mpi .NE. 0) THEN
108  WRITE(*,'("MPI_COMM_RANK returned IER =")') ierr_mpi
109  stop
110  ENDIF
111 #else
112  numprocs = 1
113  myid = 0
114 #endif
115 C
116 
117  IF(myid .eq. 0) THEN
118  WRITE(*,*) '=================================================='
119  & , '=================================================='
120  WRITE(*,*) ' This the POINCARE puncture code '
121  WRITE(*,*) ' ', trim(banner), ' Last update: ',
122  & trim(update)
123  WRITE(*,'(" Runnning on ",i5," processors")') numprocs
124  WRITE(*,*)
125  WRITE(*,*) ' Please report problems to: ',
126  & ' Don Spong (spongda@ornl.gov)'
127  & , ' or Raul Sanchez (sanchezferlr@ornl.gov).'
128  WRITE(*,*) '--------------------------------------------------'
129  & , '--------------------------------------------------'
130  ENDIF !IF(myid .eq. 0) THEN
131 
132 ! Default values for command line arguments:
133 
134  arg1 = "bfield.dat"
135  steps = 10; num_orbs = 10; scndry_steps = 50
136  nsurfs = 10; s1 = eps; s2 = 1.d0 -eps
137  angle = 0.d0
138 !
139  numargs = iargc() !for non-Cray X1 platforms; for Cray X1 use: numargs = ipxfargc()
140 
141  IF( numargs .lt. 1 .OR. trim(arg1) == '-h')THEN
142  IF(myid .eq. 0) THEN
143  WRITE(*,*)
144  print *,' Run with arguments using:'
145  WRITE(*,*)
146  WRITE(*,*) ' xpoincare file transits surfaces orbts/surf'
147  & , ' steps/period inner-surf(0->1) outer-surf(0->1) '
148  & , ' type cross-section (0->360)'
149  WRITE(*,*)
150  print *,' - FILE is the BFIELD file produced by SIESTA'
151  print *,' - transits are the number orbit transits (tor/pol)'
152  print *,' - steps/period are the steps used by integrator'
153  print *,' - inner-surf, outer-surf must be between 0 and 1.'
154  print *,' - type = poloidal or toroidal, to choose ind. var.'
155  print *,' - choose type = toroidal if BPHI nonnegative.'
156  print *,' - choose type = poloidal if BPHI =0 somewhere,'
157  print *,' and BTHETA always nonnegative.'
158  print *,' - cross-section = angle (deg) where puncture done.'
159  WRITE(*,*)
160  WRITE(*,*) ' Or run with arguments:'
161  WRITE(*,*)
162  WRITE(*,*) ' xpoincare file test_splines'
163  WRITE(*,*)
164  WRITE(*,*) ' to test splines quality with data from FILE.'
165  WRITE(*,*)
166  WRITE(*,*) '=================================================='
167  & , '=================================================='
168  ENDIF !IF(myid .eq. 0) THEN
169  stop
170  ELSE IF (numargs == 2 .AND. trim(arg2) == "test_splines") THEN
171  l_test_splines = .true.
172  ELSE
173  DO iargs = 1, numargs
174  IF (iargs == 1) CALL getarg(1, arg1)
175  IF (iargs == 2) CALL getarg(2, arg2)
176  IF (iargs == 3) CALL getarg(3, arg3)
177  IF (iargs == 4) CALL getarg(4, arg4)
178  IF (iargs == 5) CALL getarg(5, arg5)
179  IF (iargs == 6) CALL getarg(6, arg6)
180  IF (iargs == 7) CALL getarg(7, arg7)
181  IF (iargs == 8) CALL getarg(8, arg8)
182  IF (iargs == 9) CALL getarg(9, arg9)
183  END DO
184  DO iargs = numargs+1,9
185  IF (iargs == 2) arg2 = "100" !toroidal periods (50-500)
186  IF (iargs == 3) arg3 = "50" !nsurfs
187  IF (iargs == 4) arg4 = "12" !num_orbs/surf
188  IF (iargs == 5) arg5 = "40" !scndry_steps
189  IF (iargs == 6) arg6 = "0.05" !"zoom: 0.65"
190  IF (iargs == 7) arg7 = "0.95" !"zoom: 0.85"
191  IF (iargs == 8) arg8 = "toroidal"
192  IF (iargs == 9) arg9 = "0"
193  END DO
194  ENDIF
195  !for non-Cray X1 platforms. Read spline fit data, allocate arrays, open output files
196 
197  read(arg2,'(i6)') steps
198  read(arg3,'(i6)') nsurfs
199  read(arg4,'(i6)') num_orbs
200  read(arg5,'(i6)') scndry_steps
201  read(arg6,'(f8.4)') s1
202  read(arg7,'(f8.4)') s2
203  IF (s1 > s2) THEN
204  s_test = s1
205  s1 = s2
206  s2 = s1
207  ENDIF
208  IF (s1 == s2) nsurfs = 1
209  type_in = trim(arg8) ! Type = poloidal or toroidal
210  CALL tolower(type_in) ! LIBSTELL function
211  read(arg9, *) angle
212  angle = angle*(twopi/360.d0) ! Move to radians
213  angle = mod(angle, twopi) ! Make sure angle between 0. and 2*pi
214  IF (angle < 0.d0) angle = twopi + angle ! Make sure is positive
215  IF (type_in == 'poloidal') THEN
216  theta_puncture = angle
217  IF(myid .eq. 0) THEN
218  WRITE(*,'(a,1p,e10.4,a)')
219  & ' Poloidal cross section done at angle: ', angle,' rads.'
220  ENDIF !IF(myid .eq. 0) THEN
221  ELSE IF (type_in == 'toroidal') THEN
222  phi_puncture = angle
223  IF(myid .eq. 0) THEN
224  WRITE(*,'(2x,a,1p,e10.4,a)')
225  & 'Toroidal cross section done at angle: ', angle,' rads.'
226  ENDIF !IF(myid .eq. 0) THEN
227  ELSE
228  stop 'Wrong type!'
229  ENDIF
230 !
231  CALL spline_fields(arg1) ! Compute splines. Coefficients and knots in SPLINE modu
232  IF (l_test_splines) THEN ! When testing splines, do not run integration. Exit.
233  IF(myid .eq. 0) THEN
234  WRITE(*,*)
235  WRITE(*,*) ' Splines tested.'
236  WRITE(*,*)
237  CALL cpu_time(tend)
238  WRITE(*,'(a15, f12.5, a5)') ' Total time = ', tend - tstart,
239  & ' sec.'
240 
241  WRITE(*,*) '=================================================='
242  & , '=================================================='
243  ENDIF !IF(myid .eq. 0) THEN
244  stop
245 
246  ENDIF
247 
248  smin = s_knots(1) + eps ! Make sure that request is within SPLINE bounds
249  smax = (1.d0 - eps)*s_knots(is + ks)
250  thmin = u_knots(1) ! Internally, THETA runs between 0. and 2*pi
251  thmax = (1.d0 - eps)*u_knots(itht + ku)
252  phimin = v_knots(1) ! Internally, PHI runs between 0. and 2*pi
253  phimax = (1.d0 - eps) *v_knots(iphi + kv)
254  IF (s1 < smin .OR. s1 > smax) s1 = smin
255  IF (s2 < smin .OR. s2 > smax) s2 = smax
256  IF (type_in == 'toroidal') THEN
257  IF (angle < phimin .OR. angle > phimax) THEN
258  IF(myid .eq. 0) THEN
259  WRITE(*, *) 'angle requested:', angle
260  WRITE(*, *) 'PHIMIN, PHIMAX::', phimin, phimax
261  ENDIF !IF(myid .eq. 0) THEN
262  stop 'Not enough points in toroidal spline mesh!'
263  ENDIF
264  ELSEIF (type_in == 'poloidal') THEN
265  IF (angle < thmin .OR. angle > thmax) THEN
266  IF(myid .eq. 0) THEN
267  WRITE(*, *) 'angle requested:', angle
268  WRITE(*, *) 'THMIN, THMAX::', thmin, thmax
269  ENDIF !IF(myid .eq. 0) THEN
270  stop 'Not enough points in poloidal spline mesh!'
271  ENDIF
272  ENDIF
273 
274  IF(myid .eq. 0) THEN
275  WRITE(*,*)
276  WRITE(*,*) '--------------------------------------------------'
277  & , '--------------------------------------------------'
278  WRITE(*,*)
279  WRITE(*,*) ' RUN INFO:'
280  WRITE(*,*) ' ========='
281  WRITE(*,*)
282  WRITE(*,*) ' Type of integration: ', type_in
283  WRITE(*,'(" Number of toroidal periods = ",i6,2x,
284  & "using ",i6, " steps.")') steps, scndry_steps
285  WRITE(*,'(" Number of surfaces = ",i5)') nsurfs
286  WRITE(*,'(" Innermost radius of surfaces = ", f8.4)') s1
287  WRITE(*,'(" Outermost radius of surfaces = ", f8.4)') s2
288  WRITE(*,'(" Number of orbits per surface = ",i5)') num_orbs
289  ENDIF !IF(myid .eq. 0) THEN
290  local_num_orbs = nsurfs*num_orbs/numprocs
291  norbs_tot = nsurfs*num_orbs
292  num_eqns = 2*local_num_orbs; norbs = num_eqns/2 ! NUM_EQNS = Number of equations that will be advanced in time
293  IF(numprocs .gt. norbs_tot) THEN
294  WRITE(*,*)
295  WRITE(*,'("*** More processors (",i5,") than field lines ("
296  & ,i5,")***")') numprocs,norbs_tot
297  WRITE(*,*)
298  stop
299  ENDIF
300 
301  IF(numprocs*int(norbs_tot/numprocs) .ne. norbs_tot) THEN
302  WRITE(*,*)
303  WRITE(*,'("*** Number of processors (",i5,
304  & ") does not divide evenly into number of field lines ("
305  & ,i5,")***")') numprocs,norbs_tot
306  WRITE(*,*)
307  stop
308  ENDIF
309 
310  ALLOCATE(start_points(norbs_tot,2), stat=istat)
311  IF (istat .NE. 0) stop 'Allocation error'
312  ALLOCATE(theta_nb_tot(norbs_tot), stat=istat)
313  IF (istat .NE. 0) stop 'Allocation error'
314  ALLOCATE(phi_nb_tot(norbs_tot), stat=istat)
315  IF (istat .NE. 0) stop 'Allocation error'
316  ALLOCATE(s_nb_tot(norbs_tot), stat=istat)
317  IF (istat .NE. 0) stop 'Allocation error'
318  ALLOCATE(r_s_nb_tot(norbs_tot), stat=istat)
319  IF (istat .NE. 0) stop 'Allocation error'
320  ALLOCATE(z_s_nb_tot(norbs_tot), stat=istat)
321  IF (istat .NE. 0) stop 'Allocation error'
322  IF(myid .eq. 0) THEN
323  WRITE(*,'(" Total number of orbits = ", i6)') norbs
324  WRITE(*,'(" Number equations solved by LSODE = ", i6)') num_eqns
325  WRITE(*,*)
326  WRITE(*,*) '--------------------------------------------------'
327  & , '--------------------------------------------------'
328  WRITE(*,*)
329  ENDIF !IF(myid .eq. 0) THEN
330 !
331  i_truncate = 0 ! Remove .txt or .dat of input file to create output file
332  ilen = len(trim(arg1))
333  IF (arg1(ilen-3:ilen) == '.dat' .OR.
334  & arg1(ilen-3:ilen) == '.txt') i_truncate = 1
335  IF (i_truncate == 0) THEN
336  outfile = arg1(1:ilen)//'_puncture_'//trim(type_in)//'.dat'
337  ELSE
338  outfile = arg1(1:ilen-4)//'_puncture_'//trim(type_in)//'.dat'
339  ENDIF
340 
341  INQUIRE(file=trim(outfile), exist= l_ex)
342  IF (l_ex) THEN
343  OPEN(unit = 21, file = trim(outfile), status = "old"
344  & , access = "append") ! Open output file as APPEND, since then several runs can be done to complete Poincare plot
345  ELSE
346  OPEN(unit = 21, file = trim(outfile), status = "new") ! Write out headers if new file
347  IF(myid .eq. 0) THEN
348  IF (type_in == 'poloidal') THEN
349  WRITE(21,*) " & phi s r phi"
350  ELSEIF (type_in == 'toroidal') THEN
351  WRITE(21,*) " & tht s r z"
352  ENDIF
353  ENDIF !IF(myid .eq. 0) THEN
354  ENDIF
355 
356  CALL alloc_lsode ! Prepare LSODE
357 
358  icount = 0
359  DO j= 1, nsurfs ! Distribute orbit starting points for all processors
360  DO i= 1, num_orbs
361  icount = icount + 1
362  IF (nsurfs > 1) THEN
363  start_points(icount,1) = s1 + (s2-s1)*real(j-1, kind=rprec)/
364  & real(nsurfs-1, kind=rprec)
365  ELSE
366  start_points(icount,1) = s1
367  ENDIF
368  start_points(icount,2) = twopi*real(i-1, kind=rprec)/ ! If type='toroidal', even y's are THETA values; ='poloidal', they are PHI values.
369  & real(num_orbs, kind=rprec)
370  END DO
371  END DO
372 !
373  icount = 0
374 ! write(*,*) myid, myid*local_num_orbs+1,(myid+1)*local_num_orbs
375  DO i = myid*local_num_orbs+1,(myid+1)*local_num_orbs ! Distribute initial orbit locations
376  icount = icount + 1
377  y(2*icount-1) = start_points(i,1)
378  y(2*icount) = start_points(i,2)
379  ENDDO
380 !
381  IF(myid .eq. 0) THEN
382  WRITE(*,*) ' Starting integration....'
383  WRITE(*,*)
384  ENDIF !IF(myid .eq. 0) THEN
385  ang_out = 0.0d0 ! Start toroidal stepping loop
386  delta = twopi/scndry_steps ! Integration takes place with steps in toroidal angle o
387 #if defined(MPI_OPT)
388  call mpi_barrier(mpi_comm_world,ierr_mpi)
389 #endif
390  DO it = 1, steps*scndry_steps ! STEPS =No toroidal transits; SCNDRY_STEPS = steps per
391  IF(mod(it,scndry_steps) == 0) THEN
392  IF(myid .eq. 0) THEN
393  IF (type_in == 'poloidal')
394  & WRITE(*,*) ' -> Poloidal transit ',
395  & it/scndry_steps, ' completed'
396  IF (type_in == 'toroidal')
397  & WRITE(*,*) ' -> Toroidal transit ',
398  & it/scndry_steps, ' completed'
399  ENDIF !IF(myid .eq. 0) THEN
400  ENDIF
401  ang_in = ang_out ! ANG_IN = starting angle value for this loop
402  ang_out = ang_in + delta ! ANG_OUT = final angle value for this loop
403  rwork(5) = delta; iwork(6) = 5000 ! RWORK(5) = 1st step size to try; IWORK(6) = max. int
404  istate = 1; itol = 1; itask = 1; iopt = 0; jt = 10 ! ITOL = 1, same tolerance for all components; JT = 10,
405 
406  IF (type_in == 'poloidal') THEN
407  CALL lsode(rt_hand_side_dtheta, num_eqns, y, ang_in, ang_out,
408  & itol, relerr, abserr, itask, istate, iopt, rwork, lrw,
409  & iwork, liw, jacobian, jt)
410  ELSEIF (type_in == 'toroidal') THEN
411  CALL lsode(rt_hand_side_dphi, num_eqns, y, ang_in, ang_out,
412  & itol, relerr, abserr, itask, istate, iopt, rwork, lrw,
413  & iwork, liw, jacobian, jt)
414  ELSE
415  stop 'Bad type!!'
416  ENDIF
417 
418  IF (istate < 0) THEN
419  IF (istate == -1) THEN
420  WRITE(*,*) ' Excess work. Wrong JT?'
421  WRITE(*,*) ' JT = 10. Non-stiff'
422  WRITE(*,*) ' JT = 21-25. Stiff. May need JAC.'
423  ELSEIF (istate == -2) THEN
424  WRITE(*,*) ' Excess accuracy requested'
425  ELSEIF (istate == -3) THEN
426  WRITE(*,*) ' Illegal input'
427  ELSEIF (istate == -4) THEN
428  WRITE(*,*) ' Repeated error test failures'
429  ELSEIF (istate == -5) THEN
430  WRITE(*,*) ' Repeated convergence failures'
431  WRITE(*,*) ' Maybe bad Jacobian?'
432  ELSEIF (istate == -6) THEN
433  WRITE(*,*) ' Error weight became zero!'
434  ENDIF
435  WRITE(*,*) 'ANGLE done so far:', ang_in
436  ENDIF
437 !
438  IF (type_in == 'poloidal') THEN
439  DO i = 1, norbs ! Keep field lines inside spline-fit domain
440  s_test = y(2*i-1); phi_test = y(2*i)
441  CALL reset_coords(s_test, rdum, phi_test)
442  y(2*i-1) = s_test; y(2*i)= phi_test
443  ENDDO
444  ELSEIF (type_in == 'toroidal') THEN
445  DO i = 1, norbs ! Keep field lines inside spline-fit domain
446  s_test = y(2*i-1); theta_test = y(2*i)
447  CALL reset_coords(s_test, theta_test, rdum)
448  y(2*i-1) = s_test; y(2*i)= theta_test
449  ENDDO
450  ENDIF
451 
452  IF (mod(it, scndry_steps) == 0) THEN ! Write out tab-delimited puncture data
453 
454  IF (type_in == 'poloidal') THEN
455  DO i= 1, norbs
456  s_nb(i) = y(2*i-1)
457  theta_nb(i) = theta_puncture
458  phi_nb(i) = y(2*i)
459  ENDDO
460 
461  CALL dbs3vd(0, 0, 0, norbs, s_nb, theta_nb, phi_nb,
462  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi,
463  & r_coef, r_s_nb, is_uniformx, is_uniformy, is_uniformz )
464 
465 #if defined(MPI_OPT)
466  call mpi_barrier(mpi_comm_world,ierr_mpi)
467  call mpi_gather(phi_nb,norbs,mpi_real8,
468  > phi_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
469  call mpi_gather(s_nb,norbs,mpi_real8,
470  > s_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
471  call mpi_gather(r_s_nb,norbs,mpi_real8,
472  > r_s_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
473 #else
474  phi_nb_tot(:) = phi_nb(:)
475  s_nb_tot(:) = s_nb(:)
476  r_s_nb_tot(:) = r_s_nb(:)
477 #endif
478 #if defined(MPI_OPT)
479  call mpi_barrier(mpi_comm_world,ierr_mpi)
480 #endif
481  IF(myid .eq. 0) THEN
482  DO i= 1, norbs_tot
483 ! WRITE(21,'(e15.7, 3(3x, e15.7))') y(2*i), ! Pairs (phi, s) and (R, phi) at THETA=theta_puncture are written out
484 ! & y(2*i-1), r_s_nb(i), y(2*i)
485  WRITE(21,'(e15.7, 3(3x, e15.7))') phi_nb_tot(i), ! Pairs (phi, s) and (R, phi) at THETA=theta_puncture are written out
486  & s_nb_tot(i), r_s_nb_tot(i), phi_nb_tot(i)
487  END DO
488  ENDIF !IF(myid .eq. 0) THEN
489 
490  ELSEIF (type_in == 'toroidal') THEN
491 
492  DO i= 1, norbs
493  s_nb(i) = y(2*i-1)
494  theta_nb(i) = y(2*i)
495  phi_nb(i) = phi_puncture
496  ENDDO
497  CALL dbs3vd(0, 0, 0, norbs, s_nb, theta_nb, phi_nb,
498  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi,
499  & r_coef, r_s_nb, is_uniformx, is_uniformy, is_uniformz )
500 
501  CALL dbs3vd(0, 0, 0, norbs, s_nb, theta_nb, phi_nb,
502  & ks, ku, kv, s_knots, u_knots, v_knots, is, itht, iphi,
503  & z_coef, z_s_nb, is_uniformx, is_uniformy, is_uniformz )
504 
505 #if defined(MPI_OPT)
506  call mpi_barrier(mpi_comm_world,ierr_mpi)
507  call mpi_gather(theta_nb,norbs,mpi_real8,
508  > theta_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
509  call mpi_gather(s_nb,norbs,mpi_real8,
510  > s_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
511  call mpi_gather(r_s_nb,norbs,mpi_real8,
512  > r_s_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
513  call mpi_gather(z_s_nb,norbs,mpi_real8,
514  > z_s_nb_tot,norbs,mpi_real8,0,mpi_comm_world,ierr_mpi)
515 #else
516  theta_nb_tot(:) = theta_nb(:)
517  s_nb_tot(:) = s_nb(:)
518  r_s_nb_tot(:) = r_s_nb(:)
519  z_s_nb_tot(:) = z_s_nb(:)
520 #endif
521 #if defined(MPI_OPT)
522  call mpi_barrier(mpi_comm_world,ierr_mpi)
523 #endif
524  IF(myid .eq. 0) THEN
525  DO i= 1, norbs_tot
526  WRITE(21,'(e15.7, 3(3x, e15.7))') theta_nb_tot(i), ! Pairs (theta,s) and (R, Z) at PHI = phi_puncture are written out
527  & s_nb_tot(i), r_s_nb_tot(i), z_s_nb_tot(i)
528  END DO
529  ENDIF !IF(myid .eq. 0) THEN
530 
531  ENDIF
532 
533  ENDIF
534 
535  END DO
536 
537  CALL dealloc_splines ! Liberate memory
538  CALL dealloc_lsode
539 #if defined(MPI_OPT)
540  call mpi_barrier(mpi_comm_world,ierr_mpi)
541 #endif
542  DEALLOCATE(s_nb_tot,phi_nb_tot,theta_nb_tot,r_s_nb_tot,
543  & z_s_nb_tot,start_points)
544  CLOSE(21)
545 !
546  CALL cpu_time(tend)
547  IF(myid .eq. 0) THEN
548  WRITE(*,*)
549  WRITE(*,*) ' .... completed.'
550  WRITE(*,*)
551  WRITE(*,'(a15, f12.5, a5)') ' Total time = ', tend - tstart,
552  & ' sec.'
553 
554  WRITE(*,*) '=================================================='
555  & , '=================================================='
556  ENDIF !IF(myid .eq. 0) THEN
557 
558 #if defined(MPI_OPT)
559  call mpi_finalize(ierr_mpi)
560 #endif
561 
562  END PROGRAM poincare
mpi_inc
Umbrella module avoid multiple inlcudes of the mpif.h header.
Definition: mpi_inc.f:11
Generated on Thu Mar 5 2020 15:49:25 for V3FIT by   doxygen 1.8.17