V3FIT
v3fita.f
Go to the documentation of this file.
1 !*******************************************************************************
4 !
5 ! Note separating the Doxygen comment block here so detailed decription is
6 ! found in the Module not the file.
7 !
38 !*******************************************************************************
39 !*******************************************************************************
40 ! MAIN PROGRAM
41 !*******************************************************************************
42 !-------------------------------------------------------------------------------
46 !-------------------------------------------------------------------------------
47  PROGRAM v3fit
48  USE v3fit_context
49  USE profiler
50 
51  IMPLICIT NONE
52 
53 ! local variables
54  TYPE (v3fit_context_class), POINTER :: context => null()
55  REAL (rprec) :: time_start, time_end
56  INTEGER :: eq_rank
57  INTEGER :: recon_rank
58 #if defined(mpi_opt)
59  INTEGER :: error
60 
61 ! Start of executable code
62 ! The Intel version of the MPI libraries does not provide a correct value for
63 ! time until after MPI_INIT is called. Make sure this is the first think called
64 ! so that correct timing information can be used.
65  CALL mpi_init(error)
66 #endif
67 
68  CALL second0(time_start)
70 
71 ! Create a context and read in the v3fit namelist input file.
72  context => v3fit_context_construct(commandline_parser_construct())
73 
74 ! Run unit tests if the -test flag is set.
75  IF (commandline_parser_is_flag_set(context%cl_parser, &
76  & '-test')) THEN
77  my_task = 'unit_test'
78  ELSE
79  CALL v3fit_input_read_namelist( &
80  & commandline_parser_get_string(context%cl_parser, &
81  & '-file'))
82  END IF
83 
84 ! Configure parallel computation.
85  CALL config_parallelism(context)
86 
87  eq_rank = 0
88 #if defined(MPI_OPT)
89  CALL mpi_comm_rank(context%equilibrium_comm, eq_rank, error)
90 #endif
91 
92  IF (eq_rank .eq. 0) THEN
93  CALL v3fit_context_create_files(context)
94 
95  WRITE (context%recout_iou,*) ' my_task is ', my_task
96  WRITE (context%runlog_iou,*) ' my_task is ', my_task
97  END IF
98 
99  SELECT CASE (my_task)
100 
101  CASE ('equilibrium')
102  CALL task_equilibrium(context)
103 
104  CASE ('v3post', 'vmec_v3post')
105  CALL task_v3post(context)
106 
107  CASE ('reconstruct', 'reconstruct_a1')
108  CALL task_reconstruct(context)
109 
110  CASE ('child_recon')
111  CALL task_child_recon(context)
112 
113  CASE ('child_equilibrium')
114  CALL task_child_equilibrium(context)
115 
116  CASE ('unit_test')
117  CALL task_unit_tests
118 
119  CASE DEFAULT
120  CALL err_fatal('Unkown task type ' // trim(my_task))
121 
122  END SELECT
123 
124 ! Clean up memory
125  IF (eq_rank .eq. 0) THEN
126  recon_rank = 0
127 #if defined(MPI_OPT)
128  CALL mpi_comm_rank(context%reconstruction_comm, recon_rank, &
129  & error)
130 #endif
131 
132  IF (recon_rank .eq. 0) THEN
133  CALL second0(time_end)
134 
135  CALL v3fit_context_write(context)
136 
137  WRITE (context%runlog_iou,1000) time_end - time_start
138  WRITE (context%recout_iou,1000) time_end - time_start
139  WRITE (*,1000) time_end - time_start
140 
141  CALL profiler_write(context%recout_iou)
142 
143  END IF
144 
145  CALL v3fit_context_close_files(context)
146  END IF
147  CALL profiler_destruct
148 
149  CALL cleanup_parallelism(context)
150  CALL v3fit_context_destruct(context)
151 
152 #if defined(MPI_OPT)
153  CALL mpi_finalize(error)
154 #endif
155 
156 1000 FORMAT(' *** V3FIT run time (s): ',f10.2)
157 
158  END PROGRAM
159 
160 !*******************************************************************************
161 ! TASK SUBROUTINES
162 !*******************************************************************************
163 !-------------------------------------------------------------------------------
170 !-------------------------------------------------------------------------------
171  SUBROUTINE task_equilibrium(context)
173 
174  IMPLICIT NONE
175 
176 ! Declare Arguments
177  TYPE (v3fit_context_class), INTENT(inout) :: context
178 
179 ! local variables
180  INTEGER :: eq_steps
181  LOGICAL :: converged
182  INTEGER :: error
183  REAL (rprec) :: start_time
184 
185 ! Start of executable code
186  start_time = profiler_get_start_time()
187 
188  eq_steps = 1
189 
190  WRITE (*,*) ' *** In task equilibrium'
191  WRITE (context%runlog_iou,*) ' *** In task equilibrium'
192 
193 ! Initialize equilibrium
194  CALL init_equilibrium(context)
195 
196 ! Solve equilibrium
197  WRITE (*,*) ' *** Starting initial equilibrium convergence (eq)'
198  WRITE (context%runlog_iou,*) ' *** Starting initial ' // &
199  & 'equilibrium convergence'
200 
201 #if defined (MPI_OPT)
202  converged = model_converge(context%model, eq_steps, &
203  & context%runlog_iou, &
204  & context%equilibrium_comm, 'All')
205  CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0, &
206  & context%equilibrium_comm, error)
207 #else
208  converged = model_converge(context%model, eq_steps, &
209  & context%runlog_iou, 0, 'All')
210 #endif
211 
212  CALL v3fit_context_init_data(context, eq_steps)
213 
214  CALL profiler_set_stop_time('task_equilibrium', start_time)
215 
216  END SUBROUTINE
217 
218 !-------------------------------------------------------------------------------
225 !-------------------------------------------------------------------------------
226  SUBROUTINE task_v3post(context)
228 
229  IMPLICIT NONE
230 
231 ! Declare Arguments
232  TYPE (v3fit_context_class), INTENT(inout) :: context
233 
234 ! local variables
235  INTEGER :: eq_steps
236  INTEGER :: error
237  LOGICAL :: converged
238  INTEGER :: i
239  INTEGER :: last_para_signal
240  REAL (rprec) :: g2
241  REAL (rprec) :: start_time
242 
243 ! local parameters
244  REAL (rprec), DIMENSION(4), PARAMETER :: dummy_value = 0.0
245 
246 ! Start of executable code
247  start_time = profiler_get_start_time()
248 
249  eq_steps = 1
250 
251  WRITE (*,*) ' *** In task v3post'
252  WRITE (context%runlog_iou,*) ' *** In task v3post'
253 
254 ! Initialize v3post
255  CALL init_equilibrium(context)
256  CALL init_signals(context)
257  CALL init_gaussian_process(context)
258 
259 ! Solve equilibrium
260  WRITE (*,*) ' *** Starting initial equilibrium convergence (v3p)'
261  WRITE (context%runlog_iou,*) ' *** Starting initial ' // &
262  & 'equilibrium convergence'
263  converged = model_converge(context%model, eq_steps, &
264  & context%runlog_iou, &
265 #if defined(MPI_OPT)
266  & context%equilibrium_comm, 'All')
267  CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0, &
268  & context%equilibrium_comm, error)
269 #else
270  & 0, 'All')
271 #endif
272 
273 ! Set the guassian processes
274  DO i = 1, SIZE(context%gp)
275  CALL gaussp_set_profile(context%gp(i)%p, context%model)
276  END DO
277 
278 ! Calculate and display the signals.
279  WRITE (*,*) ' *** Calculating modeled signals'
280  WRITE (context%runlog_iou,*) ' *** Calculating modeled signals'
281 
282  g2 = 0.0
283 
284 ! The combination signals cannot be computed in parallel since the depend on
285 ! other signals previously computed. If combinations are used, the signals need
286 ! to be split into two loops.
287  IF (context%combination_index .eq. -1) THEN
288  last_para_signal = SIZE(context%signals)
289  ELSE
290  last_para_signal = context%combination_index - 1
291  END IF
292 
293 !$OMP PARALLEL DO
294 !$OMP& SCHEDULE(DYNAMIC)
295 !$OMP& DEFAULT(SHARED)
296 !$OMP& PRIVATE(i)
297 !$OMP& REDUCTION(+:g2)
298  DO i = 1, last_para_signal
299  g2 = g2 + signal_get_g2(context%signals(i)%p, &
300  & context%model, .false., dummy_value)
301  END DO
302 !$OMP END PARALLEL DO
303 
304  DO i = last_para_signal + 1, SIZE(context%signals)
305  g2 = g2 + context%signals(i)%p%get_g2(context%model, .false., &
306  & dummy_value)
307  END DO
308 
309  WRITE (*,1000) g2
310  WRITE (context%runlog_iou,1000) g2
311 
312  CALL v3fit_context_init_data(context, eq_steps)
313 
314  CALL profiler_set_stop_time('task_v3post', start_time)
315 
316 1000 FORMAT('g^2 = ',es12.5)
317 
318  END SUBROUTINE
319 
320 !-------------------------------------------------------------------------------
327 !-------------------------------------------------------------------------------
328  SUBROUTINE task_reconstruct(context)
330 
331  IMPLICIT NONE
332 
333 ! Declare Arguments
334  TYPE (v3fit_context_class), INTENT(inout) :: context
335 
336 ! local variables
337  INTEGER :: eq_steps
338  INTEGER :: restart_step
339  INTEGER :: i, error
340  LOGICAL :: eq_converged
341  CHARACTER (len=dir_prefix_len) :: directory
342  LOGICAL :: write_input
343  REAL (rprec) :: start_time
344 
345 ! Start of executable code
346  start_time = profiler_get_start_time()
347 
348  WRITE (*,*) ' *** In task reconstruct'
349 
350 ! Initialize the reconstruction.
351  write_input = commandline_parser_is_flag_set(context%cl_parser, &
352  & '-out')
353 
354  CALL init_equilibrium(context)
355  CALL init_signals(context)
356  CALL init_gaussian_process(context)
357  CALL init_parameters(context)
358  CALL init_reconstruction(context)
359 
360  IF (commandline_parser_is_flag_set(context%cl_parser, &
361  & '-restart')) THEN
362  WRITE (*,*) ' *** Restarting reconstruction '
363 
364  eq_steps = v3fit_context_restart(context, restart_step)
365 #if defined(MPI_OPT)
366 ! Sync the number of equilibrium steps the inital convergence took so that the'
367 ! child processes start in the same place.
368  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
369  & context%reconstruction_comm, error)
370 #endif
371  ELSE
372 ! Solve equilibrium
373  WRITE (*,*) ' *** Starting initial equilibrium ' // &
374  & 'convergence (rec)'
375  WRITE (context%runlog_iou,*) ' *** Starting initial ' // &
376  & 'equilibrium convergence'
377 
378  eq_steps = 1
379  eq_converged = reconstruction_eval_e(context%recon, &
380  & context%signals, &
381  & context%model, context%gp, &
382  & eq_steps, &
383  & context%runlog_iou, &
384 #if defined(MPI_OPT)
385  & context%equilibrium_comm)
386 #else
387  & 0)
388 #endif
389  CALL assert(eq_converged, &
390  & 'task_reconstruct no inital convergence')
391  CALL reconstruction_eval_f(context%recon, &
392  & context%derived_params, &
393  & context%model)
394 
395 ! Save the signal values from the first 0th reconstruction step.
396  DO i = 1, SIZE(context%signals)
397  context%recon%last_values(:,i) = &
398  & context%signals(i)%p%modeled
399  END DO
400 
401 #if defined(MPI_OPT)
402 ! Sync the number of equilibrium steps the inital convergence took so that the'
403 ! child processes start in the same place.
404  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
405  & context%reconstruction_comm, error)
406 #endif
407 
408  WRITE (*,*)
409  WRITE (*,1002) 0
410  WRITE (*,1003) reconstruction_get_g2(context%recon)
411 
412  WRITE (context%runlog_iou,*)
413  WRITE (context%runlog_iou,1002) 0
414  WRITE (context%runlog_iou,1003) &
415  & reconstruction_get_g2(context%recon)
416 
417  IF (write_input) THEN
418  CALL equilibrium_write_input(context%model%equilibrium, 0)
419  END IF
420 
421  CALL v3fit_context_init_data(context, eq_steps)
422 
423  restart_step = 1
424  END IF
425 
426  DO i = restart_step, SIZE(context%recon%e, 2) - 1
427  WRITE (*,*)
428  WRITE (*,1002) i
429 
430  WRITE (context%runlog_iou,*)
431  WRITE (context%runlog_iou,1002) i
432 
433  IF (reconstruction_step(context%recon, context%signals, &
434  & context%derived_params, context%locks, &
435  & context%model, context%gp, &
436  & context%params, eq_steps, &
437  & context%runlog_iou, &
438 #if defined(MPI_OPT)
439  & context%reconstruction_comm, &
440  & context%equilibrium_comm)) THEN
441 #else
442  & 0, 0)) THEN
443 #endif
444 
445 ! Step was successful. Post process the result.
446  CALL reconstruction_eval_sem(context%recon, &
447  & context%params, &
448  & context%signals, &
449  & context%derived_params)
450  CALL reconstruction_write_step(context%recon, &
451  & context%runlog_iou)
452  CALL v3fit_context_write_step_data(context, .false., &
453  & eq_steps)
454 
455  IF (write_input) THEN
456  CALL equilibrium_write_input(context%model%equilibrium, &
457  & i)
458  END IF
459 
460 ! Test for convergence.
461  IF (reconstruction_get_dg2(context%recon) .le. &
462  & context%recon_stop .or. &
463  & reconstruction_get_g2(context%recon) .eq. 0.0) THEN
464  WRITE (*,*)
465  WRITE (*,*) ' *** Equilibrium reconstructed'
466  WRITE (*,1000) context%recon_stop, &
467  & reconstruction_get_dg2(context%recon)
468 
469  WRITE (context%runlog_iou,*)
470  WRITE (context%runlog_iou,*) ' *** Equilibrium ' // &
471  & 'reconstructed'
472  WRITE (context%runlog_iou,1000) context%recon_stop, &
473  & reconstruction_get_dg2(context%recon)
474  EXIT
475  END IF
476 
477  IF (i .eq. SIZE(context%recon%e, 2) - 1) THEN
478  WRITE (*,*)
479  WRITE (*,1001) i
480 
481  WRITE (context%runlog_iou,*)
482  WRITE (context%runlog_iou,1001) i
483  END IF
484  ELSE
485  EXIT
486  END IF
487  END DO
488 
489 #if defined(MPI_OPT)
490  CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0, &
491  & context%reconstruction_comm, error)
492  CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0, &
493  & context%equilibrium_comm, error)
494 #endif
495 
496  CALL profiler_set_stop_time('task_reconstruct', start_time)
497 
498 1000 FORMAT('dg2_stop = ',es12.5,' change in g^2 = ',es12.5)
499 1001 FORMAT('nrstep 'i4,' completed. May be reconstructed.')
500 1002 FORMAT(' *** Reconstruction step ',i4)
501 1003 FORMAT('Initial g^2 = ',es12.5)
502 
503  END SUBROUTINE
504 
505 !-------------------------------------------------------------------------------
511 !-------------------------------------------------------------------------------
512  SUBROUTINE task_child_recon(context)
514 
515 ! Declare Arguments
516  TYPE (v3fit_context_class), INTENT(inout) :: context
517 
518 #if defined(MPI_OPT)
519 ! local variables
520  INTEGER :: i
521  INTEGER :: signals_created
522  INTEGER :: task
523  INTEGER :: error
524  INTEGER :: eq_steps
525  REAL (rprec) :: step_use
526  LOGICAL :: temp
527 
528 ! Start of executable code
529  CALL init_equilibrium(context)
530  CALL init_signals(context)
531  CALL init_gaussian_process(context)
532  CALL init_parameters(context)
533  CALL init_reconstruction(context)
534 
535 ! Wait for parent process to finish the initial equilibrium. Sync the number of
536 ! equilibrium steps so that the child processes match after the inital
537 ! equilibrium convergence.
538  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
539  & context%reconstruction_comm, error)
540 
541  DO
542  CALL mpi_bcast(task, 1, mpi_integer, 0, &
543  & context%reconstruction_comm, error)
544 
545  SELECT CASE (task)
546 
547  CASE (mpi_jacobian_task)
548  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
549  & context%reconstruction_comm, error)
550  CALL reconstruction_eval_jacobians(context%recon, &
551  & context%signals, context%derived_params, &
552  & context%locks, context%model, context%gp, &
553  & context%params, eq_steps, context%runlog_iou, &
554  & context%reconstruction_comm, &
555  & context%equilibrium_comm)
556 
557  CASE (mpi_sync_task)
558  CALL reconstruction_sync_svd(context%recon, &
559  & context%reconstruction_comm)
560  DO i = 1, SIZE(context%params)
561  CALL param_sync_delta(context%params(i)%p, &
562  & context%reconstruction_comm)
563  END DO
564 
565  CASE (mpi_step_task)
566  CALL mpi_bcast(step_use, 1, mpi_real8, 0, &
567  & context%reconstruction_comm, error)
568  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
569  & context%reconstruction_comm, error)
570  temp = reconstruction_try_step(context%recon, &
571  & context%signals, context%derived_params, &
572  & context%locks, context%model, context%gp, &
573  & context%params, eq_steps, step_use, &
574  & context%runlog_iou, &
575  & context%reconstruction_comm, &
576  & context%equilibrium_comm)
577 
578  CASE (mpi_quit_task)
579  CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0, &
580  & context%equilibrium_comm, error)
581  EXIT
582 
583  CASE DEFAULT
584  WRITE (*,*) 'Unknown MPI task.'
585  CALL exit(1)
586 
587  END SELECT
588  END DO
589 #endif
590 
591  END SUBROUTINE
592 
593 !-------------------------------------------------------------------------------
599 !-------------------------------------------------------------------------------
600  SUBROUTINE task_child_equilibrium(context)
602 
603 ! Declare Arguments
604  TYPE (v3fit_context_class), INTENT(inout) :: context
605 
606 #if defined(MPI_OPT)
607 ! local variables
608  INTEGER :: task
609  LOGICAL :: result
610  INTEGER :: eq_steps
611 
612 ! Start of executable code
613  eq_steps = 1
614  CALL init_equilibrium(context)
615 
616  DO
617  CALL mpi_bcast(task, 1, mpi_integer, 0, &
618  & context%equilibrium_comm, error)
619 
620  SELECT CASE (task)
621 
622  CASE (mpi_equilibrium_task)
623  CALL mpi_bcast(context%model%state_flags, 1, mpi_integer, &
624  & 0, context%equilibrium_comm, error)
625  result = equilibrium_converge(context%model%equilibrium, &
626  & eq_steps, &
627  & context%runlog_iou, &
628  & context%equilibrium_comm, &
629  & context%model%state_flags)
630 
631  CASE (mpi_mgrid_task)
632  CALL mpi_bcast(eq_steps, 1, mpi_integer, 0, &
633  & context%equilibrium_comm, error)
634  CALL equilibrium_read_vac_file(context%model%equilibrium, &
635  & eq_steps, &
636  & context%equilibrium_comm)
637 
638  CASE (mpi_quit_task)
639  EXIT
640 
641  CASE DEFAULT
642  WRITE (*,*) 'Unknown MPI task.'
643  CALL exit(1)
644 
645  END SELECT
646  END DO
647 #endif
648  END SUBROUTINE
649 
650 !-------------------------------------------------------------------------------
657 !-------------------------------------------------------------------------------
658  SUBROUTINE task_gsq_on_grid(context)
660 
661  IMPLICIT NONE
662 
663 ! Declare Arguments
664  TYPE (v3fit_context_class), INTENT(inout) :: context
665 
666 ! local variables
667  REAL (rprec) :: start_time
668 
669 ! Start of executable code
670  start_time = profiler_get_start_time()
671 
672  WRITE (*,*) ' *** In task gsq_on_grid'
673 
674 ! Initialize the reconstruction.
675  CALL init_equilibrium(context)
676  CALL init_signals(context)
677  CALL init_parameters(context)
678  CALL init_reconstruction(context)
679 
680  END SUBROUTINE
681 
682 !-------------------------------------------------------------------------------
687 !-------------------------------------------------------------------------------
688  SUBROUTINE task_unit_tests
691  USE functions
692  USE integration_path
693  USE fourier
694 
695  IMPLICIT NONE
696 
697 ! local variables
698  LOGICAL :: tests_failed
699 
700 ! Start of executable code
701  WRITE (*,*) ' *** In task unit_test'
702  tests_failed = .false.
703 
704 ! Line segment tests
705  IF (line_seg_test()) THEN
706  WRITE (*,*) 'Line Segment Test Passed'
707  ELSE
708  tests_failed = .true.
709  END IF
710 
711 ! Coordinate tests
712  IF (cood_utils_test()) THEN
713  WRITE (*,*) 'Coord Utilites Test Passed'
714  ELSE
715  tests_failed = .true.
716  END IF
717 
718 ! Profile function tests
719  IF (function_test()) THEN
720  WRITE (*,*) 'Functions Test Passed'
721  ELSE
722  tests_failed = .true.
723  END IF
724 
725 ! Integration path tests
726  IF (path_test()) THEN
727  WRITE (*,*) 'Integration Path Test Passed'
728  ELSE
729  tests_failed = .true.
730  END IF
731 
732 ! SIESTA fourier unit test.
733  IF (test_fourier()) THEN
734  WRITE (*,*) 'SIESTA Fourier Test Passed'
735  ELSE
736  tests_failed = .true.
737  END IF
738 
739 ! The CTest unit test frame work chooses a pass fail based on the exit code of
740 ! the program.
741  IF (tests_failed) THEN
742  CALL exit(1)
743  END IF
744 
745  END SUBROUTINE
746 
747 !*******************************************************************************
748 ! INITIALIZE SUBROUTINES
749 !*******************************************************************************
750 !-------------------------------------------------------------------------------
758 !-------------------------------------------------------------------------------
759  SUBROUTINE init_equilibrium(context)
761 
762  IMPLICIT NONE
763 
764 ! Declare Arguments
765  TYPE (v3fit_context_class), INTENT(inout) :: context
766 
767 ! local variables
768  LOGICAL :: force_solve
769  INTEGER :: eq_rank
770  INTEGER :: error
771  REAL (rprec) :: start_time
772 
773 ! Start of executable code
774  start_time = profiler_get_start_time()
775 
776  eq_rank = 0
777 #if defined(MPI_OPT)
778  CALL mpi_comm_rank(context%equilibrium_comm, eq_rank, error)
779 #endif
780 
781  WRITE (*,*) ' *** Constructing equilibrium model'
782  IF (eq_rank .eq. 0) THEN
783  WRITE (context%runlog_iou,*) ' *** Constructing ' // &
784  & 'equilibrium model'
785  END IF
786 
787  force_solve = &
788  & commandline_parser_is_flag_set(context%cl_parser, '-force')
789 
790  SELECT CASE (model_eq_type)
791 
792  CASE ('vmec')
793  CALL init_vmec_equilibrium(context, force_solve)
794 
795  CASE ('vacuum')
796  CALL init_vacuum_equilibrium(context, force_solve)
797 
798  CASE ('siesta')
799  CALL init_siesta_equilibrium(context, force_solve)
800 
801  CASE DEFAULT
802  WRITE (*,*) 'Equilibrium not initialized. Please ' // &
803  & 'specify a vaild equilibrium type.'
804  CALL exit(1)
805 
806  END SELECT
807 
808  CALL profiler_set_stop_time('init_equilibrium', start_time)
809 
810  END SUBROUTINE
811 
812 !-------------------------------------------------------------------------------
820 !-------------------------------------------------------------------------------
821  SUBROUTINE init_vmec_equilibrium(context, force_solve)
823  USE integration_path
824 
825  IMPLICIT NONE
826 
827 ! Declare Arguments
828  TYPE (v3fit_context_class), INTENT(inout) :: context
829  LOGICAL, INTENT(in) :: force_solve
830 
831 ! local variables
832  TYPE (pprofile_pointer), DIMENSION(:), POINTER :: sxr => null()
833  INTEGER :: i
834  LOGICAL :: not_converged
835  REAL (rprec) :: start_time
836  TYPE (emission_class), POINTER :: emission_func => null()
837  TYPE (vmec_class), POINTER :: vmec => null()
838  INTEGER :: state_flags
839 
840 ! Start of executable code
841  start_time = profiler_get_start_time()
842 
843  ALLOCATE(sxr(num_sxrem_p))
844  DO i = 1, num_sxrem_p
845  sxr(i)%p => pprofile_construct(trim(pp_sxrem_ptype_a(i)), &
846  & pp_sxrem_b_a(i,:), &
847  & pp_sxrem_as_a(i,:), &
848  & pp_sxrem_af_a(i,:))
849  END DO
850 
851 ! Create an emission object if an emission file was provided.
852  IF (trim(emission_file) .ne. '') THEN
853  emission_func => emission_construct(emission_file)
854  END IF
855 
856  state_flags = model_state_all_off
857  vmec => vmec_construct(vmec_nli_filename, vmec_wout_input, &
858  & pprofile_construct(trim(pp_ne_ptype), &
859  & pp_ne_b, pp_ne_as, &
860  & pp_ne_af), &
861  & pprofile_construct(trim(pp_te_ptype), &
862  & pp_te_b, pp_te_as, &
863  & pp_te_af), &
864  & pprofile_construct(trim(pp_ti_ptype), &
865  & pp_ti_b, pp_ti_as, &
866  & pp_ti_af), &
867  & pprofile_construct(trim(pp_ze_ptype), &
868  & pp_ze_b, pp_ze_as, &
869  & pp_ze_af), sxr, &
870  & phi_offset, z_offset, pol_rad_ratio, &
871  & context%runlog_iou, &
872 #if defined(MPI_OPT)
873  & context%equilibrium_comm, &
874  & context%reconstruction_comm, state_flags)
875 #else
876  & 0, 0, state_flags)
877 #endif
878 
879  context%model => model_construct(model_ne_type, &
880  & model_sxrem_type_a(1:num_sxrem_p), model_te_type, &
881  & model_ti_type, model_ze_type, ne_pp_unit, ne_min, te_min, &
882  & ti_min, ze_min, sxrem_min(1:num_sxrem_p), e_pressure_fraction, &
883  & emission_func, equilibrium_construct(vmec, force_solve), &
884  & sxrem_te_a, sxrem_ratio_a, ece_resonance_range, coosig_wgts, &
885  & state_flags, &
886  & sfactor_spec_fac(:minloc(sfactor_spec_imin, 1) - 1), &
887  & soffset_spec_fac(:minloc(soffset_spec_imin, 1) - 1), &
888  & path_construct(int_method, int_num_points, int_size))
889 
890  CALL profiler_set_stop_time('init_vmec_equilibrium', start_time)
891 
892  END SUBROUTINE
893 
894 !-------------------------------------------------------------------------------
902 !-------------------------------------------------------------------------------
903  SUBROUTINE init_vacuum_equilibrium(context, force_solve)
905  USE integration_path
906 
907  IMPLICIT NONE
908 
909 ! Declare Arguments
910  TYPE (v3fit_context_class), INTENT(inout) :: context
911  LOGICAL, INTENT(in) :: force_solve
912 
913 ! local variables
914  REAL (rprec) :: start_time
915  TYPE (emission_class), POINTER :: emission_func => null()
916 
917 ! Start of executable code
918  start_time = profiler_get_start_time()
919 
920 ! Create a emission object if an emission file was provided.
921  IF (trim(emission_file) .ne. '') THEN
922  emission_func => emission_construct(emission_file)
923  END IF
924 
925  context%model => model_construct(model_ne_type, &
926  & model_sxrem_type_a(1:num_sxrem_p), model_te_type, &
927  & model_ti_type, model_ze_type, ne_pp_unit, ne_min, te_min, &
928  & ti_min, ze_min, sxrem_min(1:num_sxrem_p), e_pressure_fraction, &
929  & emission_func, &
930  & equilibrium_construct(vacuum_construct(vacuum_nli_filename, &
931  & context%runlog_iou), &
932  & force_solve), &
933  & sxrem_te_a, sxrem_ratio_a, ece_resonance_range, coosig_wgts, &
934  & model_state_all_off, &
935  & sfactor_spec_fac(:minloc(sfactor_spec_imin, 1) - 1), &
936  & soffset_spec_fac(:minloc(soffset_spec_imin, 1) - 1), &
937  & path_construct(int_method, int_num_points, int_size))
938 
939  CALL profiler_set_stop_time('init_vacuum_equilibrium', start_time)
940 
941  END SUBROUTINE
942 
943 !-------------------------------------------------------------------------------
951 !-------------------------------------------------------------------------------
952  SUBROUTINE init_siesta_equilibrium(context, force_solve)
954  USE integration_path
955 
956  IMPLICIT NONE
957 
958 ! Declare Arguments
959  TYPE (v3fit_context_class), INTENT(inout) :: context
960  LOGICAL, INTENT(in) :: force_solve
961 
962 ! local variables
963  REAL (rprec) :: start_time
964  TYPE (pprofile_pointer), DIMENSION(:), POINTER :: sxr => null()
965  INTEGER :: i
966  LOGICAL :: not_converged
967  TYPE (emission_class), POINTER :: emission_func => null()
968  TYPE (vmec_class), POINTER :: vmec => null()
969  TYPE (siesta_class), POINTER :: siesta => null()
970  INTEGER :: state_flags
971 
972 ! Start of executable code
973  start_time = profiler_get_start_time()
974 
975  ALLOCATE(sxr(num_sxrem_p))
976  DO i = 1, num_sxrem_p
977  sxr(i)%p => pprofile_construct(trim(pp_sxrem_ptype_a(i)), &
978  & pp_sxrem_b_a(i,:), &
979  & pp_sxrem_as_a(i,:), &
980  & pp_sxrem_af_a(i,:))
981  END DO
982 
983 ! Create a emission object if an emission file was provided.
984  IF (trim(emission_file) .ne. '') THEN
985  emission_func => emission_construct(emission_file)
986  END IF
987 
988  state_flags = model_state_all_off
989  vmec => vmec_construct(vmec_nli_filename, vmec_wout_input, &
990  & context%runlog_iou, &
991 #if defined(MPI_OPT)
992  & context%equilibrium_comm, &
993  & context%reconstruction_comm, &
994 #else
995  & 0, 0,
996 #endif
997  & state_flags)
998 
999  siesta => siesta_construct(siesta_nli_filename, &
1000  & siesta_restart_filename, &
1001  & pprofile_construct(trim(pp_ne_ptype), &
1002  & pp_ne_b, pp_ne_as, &
1003  & pp_ne_af), &
1004  & pprofile_construct(trim(pp_te_ptype), &
1005  & pp_te_b, pp_te_as, &
1006  & pp_te_af), &
1007  & pprofile_construct(trim(pp_ti_ptype), &
1008  & pp_ti_b, pp_ti_as, &
1009  & pp_ti_af), &
1010  & sxr, phi_offset, z_offset, &
1011  & pol_rad_ratio, context%runlog_iou, &
1012 #if defined(MPI_OPT)
1013  & context%equilibrium_comm, &
1014  & context%reconstruction_comm, &
1015 #else
1016  & 0, 0, &
1017 #endif
1018  & vmec, state_flags)
1019 
1020 ! If the woutfile in the v3fit namelist input file is specifed, assume that the
1021 ! equilibrium is converged.
1022  not_converged = trim(siesta_restart_filename) .eq. ''
1023 
1024  context%model => model_construct(model_ne_type, &
1025  & model_sxrem_type_a(1:num_sxrem_p), model_te_type, &
1026  & model_te_type, model_ze_type, ne_pp_unit, ne_min, te_min, &
1027  & ti_min, ze_min, sxrem_min(1:num_sxrem_p), e_pressure_fraction, &
1028  & emission_func, equilibrium_construct_siesta(siesta, &
1029  & force_solve), &
1030  & sxrem_te_a, sxrem_ratio_a, ece_resonance_range, coosig_wgts, &
1031  & state_flags, &
1032  & sfactor_spec_fac(:minloc(sfactor_spec_imin, 1) - 1), &
1033  & soffset_spec_fac(:minloc(soffset_spec_imin, 1) - 1), &
1034  & path_construct(int_method, int_num_points, int_size))
1035 
1036  CALL profiler_set_stop_time('init_siesta_equilibrium', start_time)
1037 
1038  END SUBROUTINE
1039 
1040 !-------------------------------------------------------------------------------
1046 !-------------------------------------------------------------------------------
1047  SUBROUTINE init_signals(context)
1049 
1050  IMPLICIT NONE
1051 
1052 ! Declare Arguments
1053  TYPE (v3fit_context_class), INTENT(inout) :: context
1054 
1055 ! local variables
1056  INTEGER :: signals_created = 0
1057  INTEGER :: i
1058  INTEGER :: j
1059  REAL (rprec) :: start_time
1060 
1061 ! Start of executable code
1062  start_time = profiler_get_start_time()
1063 
1064  IF (na_s_desc .eq. 0) THEN
1065  CALL err_warn('init_signals: no signals will be created. ' // &
1066  & 'na_s_desc = 0')
1067  RETURN
1068  END IF
1069 
1070  WRITE (*,*) ' *** Constructing signals'
1071  WRITE (context%runlog_iou,*) ' *** Constructing signals'
1072 
1073  ALLOCATE(context%signals(na_s_desc))
1074 
1075 ! If a mdsig filename is given, construct the diagnostics contained inside it.
1076  IF (mdsig_list_filename .ne. '') THEN
1077  CALL init_magnetic_signals(context, signals_created)
1078  END IF
1079 
1080 ! If a sxrem filename is given, construct the diagnostics contained inside it.
1081  IF (sxrch_dot_filename .ne. '') THEN
1082  CALL init_sxrem_signals(context, signals_created)
1083  END IF
1084 
1085 ! If a intpol filename is given, construct the diagnostics contained inside it.
1086  IF (ipch_dot_filename .ne. '') THEN
1087  CALL init_intpol_signals(context, signals_created)
1088  END IF
1089 
1090 ! If a thscte filename is given, construct the diagnostics contained inside it.
1091  IF (thscte_dot_filename .ne. '') THEN
1092  CALL init_thomson_signals(context, signals_created)
1093  END IF
1094 
1095 ! If the radial extcur z position isn't negative, construct an extcurz signal.
1096  IF (extcurz_s0 .ge. 0.0) THEN
1097  CALL init_extcurz_signals(context, signals_created)
1098  END IF
1099 
1100 ! If a mse filename is given, construct the diagnostics contained inside it.
1101  IF (mse_dot_filename .ne. '') THEN
1102  CALL init_mse_signals(context, signals_created)
1103  END IF
1104 
1105 ! If an ece filename is given, construct the diagnostics contained inside it.
1106  IF (ece_dot_filename .ne. '') THEN
1107  CALL init_ece_signals(context, signals_created)
1108  END IF
1109 
1110  CALL init_limiter_signals(context, signals_created)
1111 
1112  CALL init_prior_signals(context, signals_created)
1113 
1114  IF (sxrem_ratio_dot_filename .ne. '') THEN
1115  CALL init_sxrem_ratio_signals(context, signals_created)
1116  END IF
1117 
1118  CALL init_combination_signals(context, signals_created)
1119 
1120  CALL v3fit_context_resize(context)
1121 
1122  CALL profiler_set_stop_time('init_signals', start_time)
1123 
1124  END SUBROUTINE
1125 
1126 !-------------------------------------------------------------------------------
1135 !-------------------------------------------------------------------------------
1136  SUBROUTINE init_magnetic_signals(context, signals_created)
1138  USE safe_open_mod
1139  USE v3_utilities
1140  USE ezcdf
1141  USE file_opts
1142  USE magnetic
1143 
1144  IMPLICIT NONE
1145 
1146 ! Declare Arguments
1147  TYPE (v3fit_context_class), INTENT(inout) :: context
1148  INTEGER, INTENT(inout) :: signals_created
1149 
1150 ! local variables
1151  INTEGER :: plasma_index
1152  LOGICAL :: use_point
1153  INTEGER :: mdsig_list_iou
1154  INTEGER :: mdsig_iou
1155  INTEGER :: status
1156  INTEGER :: mdsig_index
1157  CHARACTER (len=path_length) :: mdsig_path
1158  CHARACTER (len=path_length) :: mdsig_filename
1159  INTEGER :: mpi_rank
1160  class(magnetic_class), POINTER :: magnetic_object
1161  REAL (rprec) :: start_time
1162 
1163 ! Start of executable code
1164  start_time = profiler_get_start_time()
1165 
1166  status = 0
1167  mdsig_list_iou = 0
1168 
1169  mpi_rank = 0
1170 #if defined(MPI_OPT)
1171  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, status)
1172 #endif
1173  IF (mpi_rank .gt. 0) THEN
1174  IF (.not.is_absolute_path(mdsig_list_filename)) THEN
1175 ! Child jacobian processes run in a sub directory.
1176  mdsig_list_filename = build_path('..', mdsig_list_filename)
1177  END IF
1178  END IF
1179 
1180  WRITE (*,*) ' *** Constructing magnetic signals from ' // &
1181  & trim(mdsig_list_filename)
1182  WRITE (context%runlog_iou,*) ' *** Constructing magnetic ' // &
1183  & 'signals from ' // &
1184  & trim(mdsig_list_filename)
1185 
1186  CALL safe_open(mdsig_list_iou, status, trim(mdsig_list_filename), &
1187  & 'old', 'formatted')
1188  CALL assert_eq(0, status, 'init_magnetic_signals: ' // &
1189  & 'failed to open ' // trim(mdsig_list_filename))
1190 
1191 ! Get the directory of the mdsig path.
1192  mdsig_path = get_path_of_file(mdsig_list_filename)
1193 
1194 ! Need to find an index of a magnetic resonse that uses the plasma response
1195 ! function. Also need to check for the use of a point.
1196  plasma_index = 0
1197  use_point = .false.
1198 
1199  cut_comp_svd = &
1200  & commandline_parser_get_real(context%cl_parser, '-c', &
1201  & cut_comp_svd)
1202 
1203 ! Loop over each file in the list file
1204  DO
1205  READ(mdsig_list_iou, *, iostat=status) &
1206  & mdsig_index, mdsig_filename
1207 
1208 ! When the end of the file is reached, the read will return an error. Break out
1209 ! of the loop.
1210  IF (status .ne. 0) EXIT
1211 
1212 ! Check if there are two many signals.
1213  IF (signals_created + 1 .gt. v3fit_max_diagnostics) THEN
1214  CALL err_fatal('init_magnetic_signals: created signals' // &
1215  & ' exceeds v3fit_max_diagnostics')
1216  END IF
1217 
1218 ! Open the netcdf file for the mdsig.
1219  CALL cdf_open(mdsig_iou, trim(build_path(mdsig_path, &
1220  & mdsig_filename)), &
1221  & 'r', status)
1222  CALL assert_eq(0, status, 'init_magnetic_signals: ' // &
1223  & 'failed to open ', trim(mdsig_path) // &
1224  & trim(mdsig_filename))
1225 
1226 ! Construct a magnetic diagnostic signal using from the mdsig netcdf file.
1227  signals_created = signals_created + 1
1228  magnetic_object => &
1229  & magnetic_class(mdsig_iou, mag_a(signals_created), &
1230  & mag_force_coil, mag_3d_a(signals_created), &
1231  & cut_comp_svd)
1232 
1233  CALL signal_construct(magnetic_object, mdsig_iou, &
1234  & sdo_data_a(signals_created), &
1235  & sdo_sigma_a(signals_created), &
1236  & sdo_weight_a(signals_created), &
1237  & v3fit_input_find_scale_index(signals_created), &
1238  & v3fit_input_find_offset_index(signals_created))
1239 
1240  CALL cdf_close(mdsig_iou)
1241 
1242  context%signals(signals_created)%p => magnetic_object
1243 
1244  IF (plasma_index .eq. 0 .and. &
1245  & magnetic_response_use_plasma( &
1246  & magnetic_object%response)) THEN
1247  plasma_index = signals_created
1248  CALL equilibrium_set_magnetic_cache( &
1249  & context%model%equilibrium, magnetic_object%response, &
1250  & context%model%state_flags)
1251  END IF
1252 
1253  IF (.not.use_point .and. &
1254  & magnetic_response_is_point(magnetic_object%response)) THEN
1255  use_point = .true.
1256  CALL equilibrium_set_magnetic_cache( &
1257  & context%model%equilibrium, any(mag_3d_a), &
1258  & context%model%state_flags)
1259  END IF
1260 
1261 ! At lease one magnetic signal was made. Mark the index of it in the context.
1262 ! This should only be run once.
1263  IF (context%magnetic_index .eq. -1) THEN
1264  context%magnetic_index = signals_created
1265  END IF
1266  END DO
1267 
1268  CLOSE(mdsig_list_iou)
1269 
1270  CALL profiler_set_stop_time('init_magnetic_signals', start_time)
1271 
1272  END SUBROUTINE
1273 
1274 !-------------------------------------------------------------------------------
1283 !-------------------------------------------------------------------------------
1284  SUBROUTINE init_sxrem_signals(context, signals_created)
1286  USE sxrch_dot
1287  USE file_opts
1288 
1289  IMPLICIT NONE
1290 
1291 ! Declare Arguments
1292  TYPE (v3fit_context_class), INTENT(inout) :: context
1293  INTEGER, INTENT(inout) :: signals_created
1294 
1295 ! local variables
1296  INTEGER :: mpi_rank
1297  INTEGER :: error
1298  REAL (rprec) :: start_time
1299 
1300 ! Start of executable code
1301  start_time = profiler_get_start_time()
1302 
1303  mpi_rank = 0
1304 #if defined(MPI_OPT)
1305  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1306 #endif
1307  IF (mpi_rank .gt. 0) THEN
1308  IF (.not.is_absolute_path(sxrch_dot_filename)) THEN
1309 ! Child jacobian processes run in a sub directory.
1310  sxrch_dot_filename = build_path('..', sxrch_dot_filename)
1311  END IF
1312  END IF
1313 
1314  WRITE (*,*) ' *** Constructing sxrem signals from ' // &
1315  & trim(sxrch_dot_filename)
1316  WRITE (context%runlog_iou,*) ' *** Constructing sxrem ' // &
1317  & 'signals from ' // &
1318  & trim(sxrch_dot_filename)
1319 
1320 ! The setting of the sxrem_index for sxrem happens inside sxrch_dot_parse_chord
1321  CALL sxrch_dot_read(sxrch_dot_filename, context%signals, &
1322  & signals_created, sdo_data_a, sdo_sigma_a, &
1323  & sdo_weight_a, context%sxrem_index)
1324 
1325  CALL profiler_set_stop_time('init_sxrem_signals', start_time)
1326 
1327  END SUBROUTINE
1328 
1329 !-------------------------------------------------------------------------------
1338 !-------------------------------------------------------------------------------
1339  SUBROUTINE init_intpol_signals(context, signals_created)
1341  USE ipch_dot
1342  USE file_opts
1343 
1344  IMPLICIT NONE
1345 
1346 ! Declare Arguments
1347  TYPE (v3fit_context_class), INTENT(inout) :: context
1348  INTEGER, INTENT(inout) :: signals_created
1349 
1350 ! local variables
1351  LOGICAL :: use_polarimetry
1352  INTEGER :: mpi_rank
1353  INTEGER :: error
1354  REAL (rprec) :: start_time
1355 
1356 ! Start of executable code
1357  start_time = profiler_get_start_time()
1358 
1359  mpi_rank = 0
1360 #if defined(MPI_OPT)
1361  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1362 #endif
1363  IF (mpi_rank .gt. 0) THEN
1364  IF (.not.is_absolute_path(ipch_dot_filename)) THEN
1365 ! Child jacobian processes run in a sub directory.
1366  ipch_dot_filename = build_path('..', ipch_dot_filename)
1367  END IF
1368  END IF
1369 
1370  WRITE (*,*) ' *** Constructing intpol signals from ' // &
1371  & trim(ipch_dot_filename)
1372  WRITE (context%runlog_iou,*) ' *** Constructing intpol ' // &
1373  & 'signals from ' // &
1374  & trim(ipch_dot_filename)
1375 
1376 ! The setting of the intpol_index for sxrem happens inside ipch_dot_parse_chord
1377  CALL ipch_dot_read(ipch_dot_filename, context%signals, &
1378  & signals_created, sdo_data_a, sdo_sigma_a, &
1379  & sdo_weight_a, context%intpol_index, &
1380  & use_polarimetry)
1381 
1382  IF (use_polarimetry .and. &
1383  & .not.equilibrium_is_using_point( &
1384  & context%model%equilibrium)) THEN
1385  CALL equilibrium_set_magnetic_cache(context%model%equilibrium, &
1386  & .false., &
1387  & context%model%state_flags)
1388  END IF
1389 
1390  CALL profiler_set_stop_time('init_intpol_signals', start_time)
1391 
1392  END SUBROUTINE
1393 
1394 !-------------------------------------------------------------------------------
1403 !-------------------------------------------------------------------------------
1404  SUBROUTINE init_thomson_signals(context, signals_created)
1406  USE thscte_dot
1407  USE file_opts
1408 
1409  IMPLICIT NONE
1410 
1411 ! Declare Arguments
1412  TYPE (v3fit_context_class), INTENT(inout) :: context
1413  INTEGER, INTENT(inout) :: signals_created
1414 
1415 ! local variables
1416  INTEGER :: mpi_rank
1417  INTEGER :: error
1418  REAL (rprec) :: start_time
1419 
1420 ! Start of executable code
1421  start_time = profiler_get_start_time()
1422 
1423  mpi_rank = 0
1424 #if defined(MPI_OPT)
1425  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1426 #endif
1427  IF (mpi_rank .gt. 0) THEN
1428  IF (.not.is_absolute_path(thscte_dot_filename)) THEN
1429 ! Child jacobian processes run in a sub directory.
1430  thscte_dot_filename = build_path('..', thscte_dot_filename)
1431  END IF
1432  END IF
1433 
1434  WRITE (*,*) ' *** Constructing thomson signals from ' // &
1435  & trim(thscte_dot_filename)
1436  WRITE (context%runlog_iou,*) ' *** Constructing thomson ' // &
1437  & 'signals from ' // &
1438  & trim(thscte_dot_filename)
1439 
1440 ! The setting of the thomson_index for sxrem happens inside
1441 ! thscte_dot_parse_chord
1442  CALL thscte_dot_read(thscte_dot_filename, context%signals, &
1443  & signals_created, sdo_data_a, sdo_sigma_a, &
1444  & sdo_weight_a, context%thomson_index)
1445 
1446  CALL profiler_set_stop_time('init_thomson_signals', start_time)
1447 
1448  END SUBROUTINE
1449 
1450 !-------------------------------------------------------------------------------
1459 !-------------------------------------------------------------------------------
1460  SUBROUTINE init_extcurz_signals(context, signals_created)
1462  USE extcurz
1463 
1464  IMPLICIT NONE
1465 
1466 ! Declare Arguments
1467  TYPE (v3fit_context_class), INTENT(inout) :: context
1468  INTEGER, INTENT(inout) :: signals_created
1469 
1470 ! local variables
1471  INTEGER :: mpi_rank
1472  INTEGER :: error
1473  class(signal_class), POINTER :: extcurz_object
1474  REAL (rprec) :: start_time
1475 
1476 ! Start of executable code
1477  start_time = profiler_get_start_time()
1478 
1479  WRITE (*,*) ' *** Constructing extcurz signals'
1480  WRITE (context%runlog_iou,*) ' *** Constructing extcurz signals'
1481 
1482 ! Check if there are two many signals.
1483  IF (signals_created + 1 .gt. v3fit_max_diagnostics) THEN
1484  CALL err_fatal('init_extcurz_signals: created signals' // &
1485  & ' exceeds v3fit_max_diagnostics')
1486  END IF
1487 
1488  extcurz_object => extcurz_class(extcurz_s0, extcurz_u0)
1489 
1490  CALL signal_construct(extcurz_object, 'extcurz', &
1491  & 'EXTernal CURrent along Z', 'A', &
1492  & sdo_data_a(signals_created + 1), &
1493  & sdo_sigma_a(signals_created + 1), &
1494  & sdo_weight_a(signals_created + 1), &
1495  & v3fit_input_find_scale_index(signals_created + 1), &
1496  & v3fit_input_find_offset_index(signals_created + 1))
1497 
1498  context%signals(signals_created + 1)%p => extcurz_object
1499 
1500  signals_created = signals_created + 1
1501 
1502 ! At lease one extcurz signal was made. Mark the index of it in the context.
1503 ! This should only be run once.
1504  IF (context%extcurz_index .eq. -1) THEN
1505  context%extcurz_index = signals_created
1506  END IF
1507 
1508  CALL profiler_set_stop_time('init_extcurz_signals', start_time)
1509 
1510  END SUBROUTINE
1511 
1512 !-------------------------------------------------------------------------------
1521 !-------------------------------------------------------------------------------
1522  SUBROUTINE init_mse_signals(context, signals_created)
1524  USE mse_dot
1525  USE file_opts
1526 
1527  IMPLICIT NONE
1528 
1529 ! Declare Arguments
1530  TYPE (v3fit_context_class), INTENT(inout) :: context
1531  INTEGER, INTENT(inout) :: signals_created
1532 
1533 ! local variables
1534  INTEGER :: mpi_rank
1535  INTEGER :: error
1536  REAL (rprec) :: start_time
1537 
1538 ! Start of executable code
1539  start_time = profiler_get_start_time()
1540 
1541  mpi_rank = 0
1542 #if defined(MPI_OPT)
1543  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1544 #endif
1545  IF (mpi_rank .gt. 0) THEN
1546  IF (.not.is_absolute_path(mse_dot_filename)) THEN
1547 ! Child jacobian processes run in a sub directory.
1548  mse_dot_filename = build_path('..', mse_dot_filename)
1549  END IF
1550  END IF
1551 
1552  WRITE (*,*) ' *** Constructing mse signals from ' // &
1553  & trim(mse_dot_filename)
1554  WRITE (context%runlog_iou,*) ' *** Constructing mse ' // &
1555  & 'signals from ' // &
1556  & trim(mse_dot_filename)
1557 
1558 ! The setting of the mse_index for mse happens inside mse_dot_parse_chord
1559  CALL mse_dot_read(mse_dot_filename, context%signals, &
1560  & signals_created, sdo_data_a, sdo_sigma_a, &
1561  & sdo_weight_a, context%mse_index)
1562 
1563  IF (context%mse_index .gt. 0 .and. &
1564  & .not.equilibrium_is_using_point( &
1565  & context%model%equilibrium)) THEN
1566  CALL equilibrium_set_magnetic_cache(context%model%equilibrium, &
1567  & .false., &
1568  & context%model%state_flags)
1569  END IF
1570 
1571  CALL profiler_set_stop_time('init_mse_signals', start_time)
1572 
1573  END SUBROUTINE
1574 
1575 !-------------------------------------------------------------------------------
1584 !-------------------------------------------------------------------------------
1585  SUBROUTINE init_ece_signals(context, signals_created)
1587  USE ece_dot
1588  USE file_opts
1589 
1590  IMPLICIT NONE
1591 
1592 ! Declare Arguments
1593  TYPE (v3fit_context_class), INTENT(inout) :: context
1594  INTEGER, INTENT(inout) :: signals_created
1595 
1596 ! local variables
1597  INTEGER :: mpi_rank
1598  INTEGER :: error
1599  REAL (rprec) :: start_time
1600 
1601 ! Start of executable code
1602  start_time = profiler_get_start_time()
1603 
1604  mpi_rank = 0
1605 #if defined(MPI_OPT)
1606  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1607 #endif
1608  IF (mpi_rank .gt. 0) THEN
1609  IF (.not.is_absolute_path(ece_dot_filename)) THEN
1610 ! Child jacobian processes run in a sub directory.
1611  ece_dot_filename = build_path('..', ece_dot_filename)
1612  END IF
1613  END IF
1614 
1615  WRITE (*,*) ' *** Constructing ece signals from ' // &
1616  & trim(ece_dot_filename)
1617  WRITE (context%runlog_iou,*) ' *** Constructing ece ' // &
1618  & 'signals from ' // &
1619  & trim(ece_dot_filename)
1620 
1621 ! The setting of the ece_index for ece happens inside ece_dot_parse_chord
1622  CALL ece_dot_read(ece_dot_filename, context%signals, &
1623  & signals_created, sdo_data_a, sdo_sigma_a, &
1624  & sdo_weight_a, context%ece_index)
1625 
1626  IF (context%ece_index .gt. 0 .and. &
1627  & .not.equilibrium_is_using_point( &
1628  & context%model%equilibrium)) THEN
1629  CALL equilibrium_set_magnetic_cache(context%model%equilibrium, &
1630  & .false., &
1631  & context%model%state_flags)
1632  END IF
1633 
1634  CALL profiler_set_stop_time('init_ece_signals', start_time)
1635 
1636  END SUBROUTINE
1637 
1638 !-------------------------------------------------------------------------------
1646 !-------------------------------------------------------------------------------
1647  SUBROUTINE init_limiter_signals(context, signals_created)
1649  USE stel_constants, only: degree
1650  USE limiter_iso_t
1651  USE limiter_grid
1652 
1653  IMPLICIT NONE
1654 
1655 ! Declare Arguments
1656  TYPE (v3fit_context_class), INTENT(inout) :: context
1657  INTEGER, INTENT(inout) :: signals_created
1658 
1659 ! local variables
1660  INTEGER :: i
1661  CHARACTER (len=data_short_name_length) :: s_name = ''
1662  CHARACTER (len=data_name_length) :: l_name = ''
1663  INTEGER :: num_phi
1664  class(signal_class), POINTER :: limiter_obj
1665  REAL (rprec) :: start_time
1666 
1667 ! Start of executable code
1668  start_time = profiler_get_start_time()
1669 
1670  WRITE (*,*) ' *** Constructing limiter signals'
1671  WRITE (context%runlog_iou,*) ' *** Constructing limiter signals'
1672 
1673 ! For now, only limiter functions can be created.
1674  DO i = 1, n_lif
1675 ! Check if there are two many signals
1676  IF (signals_created + 1 .gt. v3fit_max_signals) THEN
1677  CALL err_fatal('init_limiter_signals: created signals' // &
1678  & ' exceeds v3fit_max_signals')
1679  END IF
1680 
1681 ! Check if there are two many limiter signals.
1682  IF (i .gt. v3fit_max_limiters) THEN
1683  CALL err_fatal('init_limiter_signals: created signals' // &
1684  & ' exceeds v3fit_max_limits')
1685  END IF
1686 
1687  s_name = trim(signal_make_short_name('edge_lim_', i))
1688  WRITE (l_name,'(a,i4)') 'edge limiter iso function ', i
1689 
1690  num_phi = max(min(v3fit_max_lif_size, n_phi_lif(i)), 1)
1691 
1692 ! Check if the sdo_sigma is zero and replace it lif_sigma instead.
1693  IF (sdo_sigma_a(signals_created + 1) .eq. 0.0) THEN
1694  sdo_sigma_a(signals_created + 1) = lif_sigma(i)
1695  END IF
1696 
1697  limiter_obj => &
1698  & limiter_iso_class(lif_arz(i,:,:), lif_rc(i), lif_zc(i), 20, &
1699  & lif_phi_degree(i,1:num_phi)*degree, &
1700  & lif_on_edge(i))
1701 
1702  CALL signal_construct(limiter_obj, s_name, l_name, 'm', &
1703  & sdo_data_a(signals_created + 1), &
1704  & sdo_sigma_a(signals_created + 1), &
1705  & sdo_weight_a(signals_created + 1), &
1706  & v3fit_input_find_scale_index(signals_created + 1), &
1707  & v3fit_input_find_offset_index(signals_created + 1))
1708 
1709  context%signals(signals_created + 1)%p => limiter_obj
1710 
1711  signals_created = signals_created + 1
1712 
1713 ! At lease one limiter signal was made. Mark the index of it in the context.
1714 ! This should only be run once.
1715  IF (context%limiter_index .eq. -1) THEN
1716  context%limiter_index = signals_created
1717  END IF
1718 
1719  END DO
1720 
1721 ! Add grid limiters to the end.
1722  IF (limiter_grid_file .ne. '') THEN
1723  s_name = trim(signal_make_short_name('edge_lim_', i))
1724  WRITE (l_name,'(a,i4)') 'edge limiter polygon ', i
1725 
1726  limiter_obj => limiter_grid_class(limiter_grid_file, &
1727  & lif_on_edge(n_lif + 1))
1728 
1729  CALL signal_construct( limiter_obj, s_name, l_name, 'm', &
1730  & sdo_data_a(signals_created + 1), &
1731  & sdo_sigma_a(signals_created + 1), &
1732  & sdo_weight_a(signals_created + 1), &
1733  & v3fit_input_find_scale_index(signals_created + 1), &
1734  & v3fit_input_find_offset_index(signals_created + 1))
1735 
1736  context%signals(signals_created + 1)%p => limiter_obj
1737 
1738  signals_created = signals_created + 1
1739 
1740 ! At lease one limiter signal was made. Mark the index of it in the context.
1741 ! This should only be run once.
1742  IF (context%limiter_index .eq. -1) THEN
1743  context%limiter_index = signals_created
1744  END IF
1745 
1746  END IF
1747 
1748  CALL profiler_set_stop_time('init_limiter_signals', start_time)
1749 
1750  END SUBROUTINE
1751 
1752 !-------------------------------------------------------------------------------
1760 !-------------------------------------------------------------------------------
1761  SUBROUTINE init_prior_signals(context, signals_created)
1763  USE prior_gaussian
1764 
1765  IMPLICIT NONE
1766 
1767 ! Declare Arguments
1768  TYPE (v3fit_context_class), INTENT(inout) :: context
1769  INTEGER, INTENT(inout) :: signals_created
1770 
1771 ! local variables
1772  class(signal_class), POINTER :: prior_obj
1773  INTEGER :: i
1774  REAL (rprec) :: start_time
1775 
1776 ! Start of executable code
1777  start_time = profiler_get_start_time()
1778 
1779  WRITE (*,*) ' *** Constructing prior signals'
1780  WRITE (context%runlog_iou,*) ' *** Constructing prior signals'
1781 
1782 ! For now, only limiter functions can be created.
1783  DO i = 1, n_prior
1784 ! Check if there are two many signals
1785  IF (signals_created + 1 .gt. v3fit_max_signals) THEN
1786  CALL err_fatal('init_prior_signals: created signals' // &
1787  & ' exceeds v3fit_max_signals')
1788  END IF
1789 
1790 ! Check if there are two many prior signals.
1791  IF (i .gt. v3fit_max_priors) THEN
1792  CALL err_fatal('init_prior_signals: created signals' // &
1793  & ' exceeds v3fit_max_priors')
1794  END IF
1795 
1796  prior_obj => prior_gaussian_class(context%model, &
1797  & prior_param_name(i), &
1798  & prior_indices(i,:))
1799 
1800  CALL signal_construct(prior_obj, &
1801  & trim(signal_make_short_name('prior_', i)), &
1802  & prior_name(i), prior_units(i), &
1803  & sdo_data_a(signals_created + 1), &
1804  & sdo_sigma_a(signals_created + 1), &
1805  & sdo_weight_a(signals_created + 1), &
1806  & v3fit_input_find_scale_index(signals_created + 1), &
1807  & v3fit_input_find_offset_index(signals_created + 1))
1808 
1809  context%signals(signals_created + 1)%p => prior_obj
1810 
1811  signals_created = signals_created + 1
1812 
1813 ! At lease one prior signal was made. Mark the index of it in the context.
1814 ! This should only be run once.
1815  IF (context%prior_gaussian_index .eq. -1) THEN
1816  context%prior_gaussian_index = signals_created
1817  END IF
1818 
1819  END DO
1820 
1821  CALL profiler_set_stop_time('init_prior_signals', start_time)
1822 
1823  END SUBROUTINE
1824 
1825 !-------------------------------------------------------------------------------
1834 !-------------------------------------------------------------------------------
1835  SUBROUTINE init_sxrem_ratio_signals(context, signals_created)
1837  USE sxrem_ratio_dot
1838  USE file_opts
1839 
1840  IMPLICIT NONE
1841 
1842 ! Declare Arguments
1843  TYPE (v3fit_context_class), INTENT(inout) :: context
1844  INTEGER, INTENT(inout) :: signals_created
1845 
1846 ! local variables
1847  INTEGER :: mpi_rank
1848  INTEGER :: error
1849  REAL (rprec) :: start_time
1850 
1851 ! Start of executable code
1852  start_time = profiler_get_start_time()
1853 
1854  mpi_rank = 0
1855 #if defined(MPI_OPT)
1856  CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1857 #endif
1858  IF (mpi_rank .gt. 0) THEN
1859  IF (.not.is_absolute_path(sxrem_ratio_dot_filename)) THEN
1860 ! Child jacobian processes run in a sub directory.
1861  sxrem_ratio_dot_filename = &
1862  & build_path('..', sxrem_ratio_dot_filename)
1863  END IF
1864  END IF
1865 
1866  WRITE (*,*) ' *** Constructing sxrem ratio signals from ' // &
1867  & trim(sxrem_ratio_dot_filename)
1868  WRITE (context%runlog_iou,*) ' *** Constructing intpol ' // &
1869  & 'signals from ' // &
1870  & trim(sxrem_ratio_dot_filename)
1871 
1872 ! The setting of the intpol_index for sxrem happens inside ipch_dot_parse_chord
1873  CALL sxrem_ratio_dot_read(sxrem_ratio_dot_filename, &
1874  & context%signals, signals_created, &
1875  & sdo_data_a, sdo_sigma_a, sdo_weight_a, &
1876  & context%sxrem_ratio_index)
1877 
1878  CALL profiler_set_stop_time('init_sxrem_ratio_signals', &
1879  & start_time)
1880 
1881  END SUBROUTINE
1882 
1883 !-------------------------------------------------------------------------------
1889 !-------------------------------------------------------------------------------
1890  SUBROUTINE init_combination_signals(context, signals_created)
1892 
1893  IMPLICIT NONE
1894 
1895 ! Declare Arguments
1896  TYPE (v3fit_context_class), INTENT(inout) :: context
1897  INTEGER, INTENT(inout) :: signals_created
1898 
1899 ! local variables
1900  INTEGER :: i
1901  INTEGER :: j
1902  class(combination_class), POINTER :: combination_signal
1903  REAL (rprec) :: start_time
1904 
1905 ! Start of executable code
1906  start_time = profiler_get_start_time()
1907 
1908  WRITE (*,*) ' *** Constructing combination signals'
1909  WRITE (context%runlog_iou,*) ' *** Constructing ' // &
1910  & 'combination signals'
1911 
1912  DO i = 1, n_coosig
1913 ! Check if there are two many signals
1914  IF (signals_created + 1 .gt. v3fit_max_signals) THEN
1915  CALL err_fatal('init_combination_signals: created ' // &
1916  & 'signals exceeds v3fit_max_signals')
1917  END IF
1918 
1919 ! Check if there are two many combination signals.
1920  IF (i .gt. v3fit_max_combinations) THEN
1921  CALL err_fatal('init_combination_signals: created ' // &
1922  & 'signals exceeds v3fit_max_combinations')
1923  END IF
1924 
1925 ! Search for the last index.
1926  combination_signal => &
1927  & combination_construct(n_sig_coosig(i), coosig_type(i), &
1928  & coosig_wgts_id(i))
1929  DO j = 1, n_sig_coosig(i)
1930  IF (coosig_indices(i,j) .gt. signals_created) THEN
1931  CALL err_fatal('init_combination_signals: cannot ' // &
1932  & 'create combination signal with ' // &
1933  & 'uninitalized signal')
1934  END IF
1935 
1936  CALL combination_set_signal(combination_signal, &
1937  & context%signals(coosig_indices(i,j))%p, &
1938  & coosig_coeff(i,j), j)
1939  END DO
1940 
1941  CALL signal_construct(combination_signal, coosig_name(i), &
1942  & coosig_name(i), coosig_units(i), &
1943  & sdo_data_a(signals_created + 1), &
1944  & sdo_sigma_a(signals_created + 1), &
1945  & sdo_weight_a(signals_created + 1), &
1946  & v3fit_input_find_scale_index(signals_created + 1), &
1947  & v3fit_input_find_offset_index(signals_created + 1))
1948 
1949  context%signals(signals_created + 1)%p => combination_signal
1950 
1951  signals_created = signals_created + 1
1952 
1953 ! At lease one combination signal was made. Mark the index of it in the
1954 ! context. This should only be run once.
1955  IF (context%combination_index .eq. -1) THEN
1956  context%combination_index = signals_created
1957  END IF
1958 
1959  END DO
1960 
1961  CALL profiler_set_stop_time('init_combination_signals', &
1962  & start_time)
1963 
1964  END SUBROUTINE
1965 
1966 !-------------------------------------------------------------------------------
1973 !-------------------------------------------------------------------------------
1974  SUBROUTINE init_parameters(context)
1976 
1977  IMPLICIT NONE
1978 
1979 ! Declare Arguments
1980  TYPE (v3fit_context_class), INTENT(inout) :: context
1981 
1982 ! local variables
1983  INTEGER :: i, j
1984  INTEGER, DIMENSION(:,:), ALLOCATABLE :: indices
1985  REAL (rprec) :: start_time
1986 
1987 ! Start of executable code
1988  start_time = profiler_get_start_time()
1989 
1990  WRITE (*,*) ' *** Constructing reconstruction parameters'
1991  WRITE (context%runlog_iou,*) ' *** Constructing ' // &
1992  & 'reconstruction parameters'
1993 
1994 ! Initialize the reconstruction parameters.
1995  ALLOCATE(context%params(n_rp))
1996  DO i = 1, n_rp
1997  context%params(i)%p => param_construct(context%model, &
1998  & rp_type(i), &
1999  & (/ rp_index(i), &
2000  & rp_index2(i) /), &
2001  & rp_vrnc(i), &
2002  & rp_range_type(i,:), &
2003  & rp_range_index(i,:,:), &
2004  & rp_range_value(i,:), &
2005  & SIZE(context%signals), &
2006  & n_rp)
2007  END DO
2008 
2009 ! Initialize the derived parameters.
2010  ALLOCATE(context%derived_params(n_dp))
2011  DO i = 1, n_dp
2012  context%derived_params(i)%p => param_construct(context%model, &
2013  & dp_type(i), &
2014  & dp_index(i,:), &
2015  & n_dp)
2016 
2017  END DO
2018 
2019 ! Initialize the locking parameters.
2020  ALLOCATE(context%locks(n_lp))
2021 
2022  DO i = 1, n_lp
2023 ! Find the number of parameters to lock to.
2024  DO j = 1, SIZE(lp_sets,2)
2025  IF (lp_sets(i,j) .eq. '') THEN
2026  EXIT
2027  END IF
2028  END DO
2029  j = min(j - 1, SIZE(lp_sets,2))
2030 
2031  CALL assert(j .ne. 0, 'Locking parameters require at ' // &
2032  & 'least one set.')
2033 
2034  ALLOCATE(indices(j, data_max_indices))
2035  indices(:,1) = lp_sets_index(i,1:j)
2036  indices(:,2) = lp_sets_index2(i,1:j)
2037 
2038  context%locks(i)%p => &
2039  & param_construct(context%model, lp_type(i), &
2040  & (/ lp_index(i), lp_index2(i) /), &
2041  & lp_sets(i,1:j), indices, &
2042  & lp_sets_coeff(i,1:j), &
2043  & context%equilibrium_comm)
2044 
2045  DEALLOCATE(indices)
2046  END DO
2047 
2048  CALL profiler_set_stop_time('init_parameters', start_time)
2049 
2050  END SUBROUTINE
2051 
2052 !-------------------------------------------------------------------------------
2059 !-------------------------------------------------------------------------------
2060  SUBROUTINE init_gaussian_process(context)
2062 
2063  IMPLICIT NONE
2064 
2065 ! Declare Arguments
2066  TYPE (v3fit_context_class), INTENT(inout) :: context
2067 
2068 ! local variables
2069  INTEGER :: i
2070  INTEGER :: j
2071  REAL (rprec) :: start_time
2072 
2073 ! Start of executable code
2074  start_time = profiler_get_start_time()
2075 
2076  ALLOCATE(context%gp(n_gp))
2077 
2078  DO i = 1, n_gp
2079  context%gp(i)%p => gaussp_construct(context%model, &
2080  & n_gp_signal(i), &
2081  & gp_model_type(i), &
2082  & gp_model_index(i), &
2083  & gp_param_vrnc(i,:), &
2084  & gp_tolerance(i), &
2085  & gp_cholesky_fact(i))
2086 
2087  DO j = 1, n_gp_signal(i)
2088  CALL gaussp_set_signal(context%gp(i)%p, &
2089  & context%signals(gp_signal_indices(i,j))%p, j)
2090  END DO
2091  END DO
2092 
2093  CALL profiler_set_stop_time('init_gaussian_process', start_time)
2094 
2095  END SUBROUTINE
2096 
2097 !-------------------------------------------------------------------------------
2101 !-------------------------------------------------------------------------------
2102  SUBROUTINE init_reconstruction(context)
2104 
2105  IMPLICIT NONE
2106 
2107 ! Declare Arguments
2108  TYPE (v3fit_context_class), INTENT(inout) :: context
2109 
2110 ! local variables
2111  INTEGER :: last_para_signal
2112  REAL (rprec) :: start_time
2113 
2114 ! Start of executable code
2115  start_time = profiler_get_start_time()
2116 
2117 ! The combination signals cannot be computed in parallel since the depend on
2118 ! other signals previously computed. If combinations are used, the signals need
2119 ! to be split into two loops.
2120  IF (context%combination_index .eq. -1) THEN
2121  last_para_signal = SIZE(context%signals)
2122  ELSE
2123  last_para_signal = context%combination_index - 1
2124  END IF
2125 
2126 ! Override from command line if set.
2127  use_central_diff = use_central_diff .or. &
2128  & commandline_parser_is_flag_set(context%cl_parser, '-c_diff')
2129 
2130  context%recon_stop = dg2_stop
2131  context%recon => reconstruction_construct(nrstep, &
2132  & SIZE(context%signals), &
2133  & SIZE(context%derived_params), &
2134  & SIZE(context%params), &
2135  & step_type, astep_max, cut_svd, cut_eff, &
2136  & cut_marg_eff, cut_delta_a, cut_dg2, &
2137  & last_para_signal, cut_inv_svd, &
2138  & use_central_diff)
2139 
2140  CALL profiler_set_stop_time('init_reconstruction', start_time)
2141 
2142  END SUBROUTINE
2143 
2144 !-------------------------------------------------------------------------------
2168 !-------------------------------------------------------------------------------
2169 ! When using openmp the !$ allows conditional computation based on if openmp is
2170 ! used or not.
2171  SUBROUTINE config_parallelism(context)
2173  USE file_opts
2174 !$ USE omp_lib
2175 #if defined(MPI_OPT)
2176  USE data_parameters
2177 #endif
2178 
2179  IMPLICIT NONE
2180 
2181 ! Declare Arguments
2182  TYPE (v3fit_context_class), INTENT(inout) :: context
2183 
2184 ! local variables
2185  INTEGER :: num_threads
2186 #if defined(MPI_OPT)
2187  INTEGER :: error
2188  INTEGER :: num_recon_processes
2189  INTEGER :: num_eq_processes
2190  INTEGER :: i
2191  INTEGER :: mpi_rank
2192  INTEGER :: mpi_size
2193  INTEGER :: eq_rank
2194  INTEGER :: eq_size
2195  INTEGER :: recon_rank
2196  CHARACTER (len=dir_prefix_len) :: directory_name
2197  CHARACTER (len=path_length) :: temp_string
2198 #endif
2199  REAL (rprec) :: start_time
2200 
2201 ! Start of executable code
2202  start_time = profiler_get_start_time()
2203 
2204  WRITE (*,*) ' *** Configuring parallelism.'
2205 
2206 !-- OpenMP ---------------------------------------------------------------------
2207 
2208 ! Configure the number of threeds to use.
2209 !$ num_threads = commandline_parser_get_integer(context%cl_parser, &
2210 !$ & '-para', 1)
2211 
2212 !$ IF (num_threads .gt. 0) THEN
2213 !$ CALL OMP_SET_NUM_THREADS(num_threads)
2214 !$ END IF
2215 !$OMP PARALLEL
2216 !$ num_threads = OMP_GET_MAX_THREADS()
2217 !$OMP END PARALLEL
2218 !$ WRITE (*,1000) num_threads
2219 
2220 !-- MPI ------------------------------------------------------------------------
2221 
2222 #if defined (MPI_OPT)
2223 ! Configure num processes to use for MPI
2224  CALL mpi_comm_size(mpi_comm_world, mpi_size, error)
2225  CALL mpi_comm_rank(mpi_comm_world, mpi_rank, error)
2226 
2227 ! Initialize the intra communicator to the world communicator. This is done so
2228 ! global_comm has is value even when num_processes is 1 or the processes were
2229 ! created externally.
2230  context%global_comm = mpi_comm_world
2231 
2232 ! MPI parallelization is only implemented for specific tasks. If the root task
2233 ! is not a reconstruction task stop the child processes.
2234  SELECT CASE (my_task)
2235 
2236  CASE ('reconstruct', 'reconstruct_a1')
2237  IF (commandline_parser_is_flag_set(context%cl_parser, &
2238  & '-serial')) THEN
2239  num_recon_processes = 1
2240  num_eq_processes = mpi_size
2241  ELSE IF (mpi_size .le. n_rp) THEN
2242  num_recon_processes = mpi_size
2243  num_eq_processes = 1
2244  ELSE
2245  num_recon_processes = n_rp
2246  num_eq_processes = mpi_size/n_rp &
2247  & + mod(mpi_size, n_rp)
2248  END IF
2249 
2250  CASE ('v3post', 'vmec_v3post', 'equilibrium')
2251  num_recon_processes = 1
2252  num_eq_processes = mpi_size
2253 
2254  CASE DEFAULT
2255  num_recon_processes = 1
2256  IF (mpi_rank .gt. 0) THEN
2257  stop 'MPI is only supported for the reconstruct task.'
2258  END IF
2259 
2260  END SELECT
2261 
2262 ! Create sub working directories for the processes. When running in default
2263 ! mode one process can use the main directory. Otherwise create one directory
2264 ! for each sub process and spawn a child process in it. If this point is
2265 ! reached and the size is greater than 1, we can assume that default was used.
2266  temp_string = commandline_parser_get_string(context%cl_parser, &
2267  & '-file')
2268 
2269  IF (mpi_rank .eq. 0) THEN
2270  DO i = 2, num_recon_processes
2271 ! Create the directory and copy the namelist input files over.
2272  CALL create_directory(process_dir(i), error)
2273 
2274  IF (.not.is_absolute_path(temp_string)) THEN
2275  CALL copy_file(temp_string, &
2276  & build_path(process_dir(i), temp_string), &
2277  & error)
2278  END IF
2279 
2280 ! Copy the equilibrium inputs.
2281  IF (.not.is_absolute_path(vmec_nli_filename)) THEN
2282  CALL copy_file(vmec_nli_filename, &
2283  & build_path(process_dir(i), &
2284  & vmec_nli_filename), error)
2285  END IF
2286  IF (.not.is_absolute_path(vmec_wout_input)) THEN
2287  CALL copy_file(vmec_wout_input, &
2288  & build_path(process_dir(i), &
2289  & vmec_wout_input), error)
2290  END IF
2291  IF (.not.is_absolute_path(vacuum_nli_filename)) THEN
2292  CALL copy_file(vacuum_nli_filename, &
2293  & build_path(process_dir(i), &
2294  & vacuum_nli_filename), error)
2295  END IF
2296  IF (.not.is_absolute_path(siesta_nli_filename)) THEN
2297  CALL copy_file(siesta_nli_filename, &
2298  & build_path(process_dir(i), &
2299  & siesta_nli_filename), error)
2300  END IF
2301  IF (.not.is_absolute_path(siesta_restart_filename)) THEN
2302  CALL copy_file(siesta_restart_filename, &
2303  & build_path(process_dir(i), &
2304  & siesta_restart_filename), &
2305  & error)
2306  END IF
2307  END DO
2308  END IF
2309 
2310 ! Wait for all processes to reach this point.
2311  CALL mpi_barrier(context%global_comm, error)
2312 
2313 ! Partition the equilibrium communicators.
2314  recon_rank = mod(mpi_rank, max(1, num_recon_processes))
2315  CALL mpi_comm_split(context%global_comm, recon_rank, 0, &
2316  & context%equilibrium_comm, error)
2317  CALL mpi_comm_rank(context%equilibrium_comm, eq_rank, error)
2318 
2319 ! Take the root process of each equilibrium communicator and partition them
2320 ! from global to the reconstuction communicator.
2321  IF (eq_rank .gt. 0) THEN
2322  CALL mpi_comm_split(context%global_comm, mpi_undefined, 0, &
2323  & context%reconstruction_comm, error)
2324  ELSE
2325  CALL mpi_comm_split(context%global_comm, 0, recon_rank, &
2326  & context%reconstruction_comm, error)
2327  END IF
2328 
2329 ! Change the current directory to the reconstruction directory.
2330  IF (recon_rank .gt. 0) THEN
2331  CALL change_directory(process_dir(recon_rank + 1), error)
2332  END IF
2333 
2334  IF (mpi_rank .gt. 0) THEN
2335  IF (eq_rank .gt. 0) THEN
2336  my_task = 'child_equilibrium'
2337  ELSE
2338  my_task = 'child_recon'
2339  END IF
2340  END IF
2341 
2342  WRITE (*,1001) mpi_rank + 1, recon_rank + 1, eq_rank + 1, &
2343  & trim(my_task)
2344 
2345  CALL mpi_barrier(context%global_comm, error)
2346 #endif
2347 
2348  CALL profiler_set_stop_time('config_parallelism', start_time)
2349 
2350 1000 FORMAT('Using ',i4,' threads.')
2351 1001 FORMAT('Configured global process ',i5,', recon process ',i5, &
2352  & ', eq process ',i5,' task ',a)
2353 
2354  END SUBROUTINE
2355 
2356 !-------------------------------------------------------------------------------
2362 !-------------------------------------------------------------------------------
2363  SUBROUTINE cleanup_parallelism(context)
2365  USE file_opts
2366 #if defined(MPI_OPT)
2367  USE data_parameters
2368 #endif
2369 
2370  IMPLICIT NONE
2371 
2372 ! Declare Arguments
2373  TYPE (v3fit_context_class), INTENT(inout) :: context
2374 
2375 ! local variables
2376 #if defined(MPI_OPT)
2377  INTEGER :: error
2378  INTEGER :: mpi_rank
2379  INTEGER :: mpi_size
2380  INTEGER :: recon_size
2381  INTEGER :: i
2382 
2383 ! Start of executable code
2384  IF (context%reconstruction_comm .eq. mpi_comm_null) THEN
2385  RETURN
2386  END IF
2387 
2388  CALL mpi_comm_rank(context%global_comm, mpi_rank, error)
2389 
2390  CALL mpi_comm_size(context%reconstruction_comm, recon_size, &
2391  & error)
2392  CALL mpi_barrier(context%reconstruction_comm, error)
2393 
2394 
2395  IF (mpi_rank .eq. 0) THEN
2396 ! Clean up sub working directories for the processes. Parent directory should
2397 ! be ignored.
2398  DO i = 2, recon_size
2399  CALL delete_directory(process_dir(i), error)
2400  END DO
2401  END IF
2402 #endif
2403 
2404  END SUBROUTINE
init_magnetic_signals
subroutine init_magnetic_signals(context, signals_created)
Constructs magnetic signals.
Definition: v3fita.f:1137
v3fit_context::v3fit_context_close_files
subroutine v3fit_context_close_files(this)
Close output files.
Definition: v3fit_context.f:419
extcurz
Represents a signal obtained by integrating around the magnetic field to get the current enclosed in ...
Definition: extcurz.f:15
task_unit_tests
subroutine task_unit_tests
Run unit tests.
Definition: v3fita.f:689
file_opts::change_directory
subroutine change_directory(path, error)
Change working directory.
Definition: file_opts.f:414
sxrch_dot
Module for opening and reading a 'sxrch.' file. The file format for these files are documented in Sof...
Definition: sxrch_dot.f:108
init_equilibrium
subroutine init_equilibrium(context)
Initialize an equilibrium.
Definition: v3fita.f:760
profiler
Defines functions for measuring an tabulating performance of function and subroutine calls....
Definition: profiler.f:13
coordinate_utilities
Module is part of the LIBSTELL. This modules containes code to convert from different coordinate syst...
Definition: coordinate_utilities.f:12
init_extcurz_signals
subroutine init_extcurz_signals(context, signals_created)
Constructs the extcurz signal.
Definition: v3fita.f:1461
profiler::profiler_destruct
subroutine profiler_destruct()
Deconstruct a profiler.
Definition: profiler.f:91
thscte_dot::thscte_dot_read
subroutine thscte_dot_read(thscte_file, signals, signals_created, observed, sigma, weight, first_index)
Read an thomson scattering diagnostic dot file.
Definition: thscte_dot.f:103
file_opts::copy_file
subroutine copy_file(file_source, file_dest, error)
Copies the source file to the destination.
Definition: file_opts.f:78
integration_path::path_construct
Construction interface using either path_construct_int or path_construct_vertex.
Definition: integration_path.f:75
mse_dot
Module for opening and reading a 'mse.' file. The file format for these files are documented in Motio...
Definition: mse_dot.f:129
task_child_recon
subroutine task_child_recon(context)
Child task to compute jacobian.
Definition: v3fita.f:513
config_parallelism
subroutine config_parallelism(context)
Configure Parallelism.
Definition: v3fita.f:2172
magnetic
Defines the base class of the type magnetic_class.
Definition: magnetic.f:13
v3fit_context
Defines a v3fit_context_class object to contain all the memory for running v3fit. Contains methods to...
Definition: v3fit_context.f:12
v3_utilities::assert_eq
Definition: v3_utilities.f:62
init_siesta_equilibrium
subroutine init_siesta_equilibrium(context, force_solve)
Initialize a siests equilibrium.
Definition: v3fita.f:953
mse_dot::mse_dot_read
subroutine mse_dot_read(mse_file, signals, signals_created, observed, sigma, weight, first_index)
Read an motional stark effect diagnostic dot file.
Definition: mse_dot.f:159
v3fit_context::v3fit_context_restart
integer function v3fit_context_restart(this, current_step)
Restart the reconstruction.
Definition: v3fit_context.f:1285
file_opts::get_path_of_file
character(len=path_length) function get_path_of_file(file)
Returns the directory of a file.
Definition: file_opts.f:354
prior_gaussian::prior_gaussian_class
Base class representing a prior_guassian signal.
Definition: prior_gaussian.f:38
magnetic::magnetic_class
Base class representing a magnetic signal.
Definition: magnetic.f:39
init_sxrem_ratio_signals
subroutine init_sxrem_ratio_signals(context, signals_created)
Constructs sxrem_ratio signals.
Definition: v3fita.f:1836
sxrem_ratio_dot
Module for opening and reading a 'sxrem_ratio.' file. The file format for these files are documented ...
Definition: sxrem_ratio_dot.f:73
sxrch_dot::sxrch_dot_read
subroutine sxrch_dot_read(sxrch_file, signals, signals_created, observed, sigma, weight, first_index)
Read a soft x-ray diagnostic dot file.
Definition: sxrch_dot.f:139
init_mse_signals
subroutine init_mse_signals(context, signals_created)
Constructs mse signals.
Definition: v3fita.f:1523
v3fit_context::v3fit_context_create_files
subroutine v3fit_context_create_files(this)
Create output files.
Definition: v3fit_context.f:339
task_gsq_on_grid
subroutine task_gsq_on_grid(context)
Computes the g^2 value on a grid.
Definition: v3fita.f:659
init_ece_signals
subroutine init_ece_signals(context, signals_created)
Constructs ece signals.
Definition: v3fita.f:1586
fourier
Module contains subroutines for computing FFT on parallel radial intervals. Converts quantities betwe...
Definition: fourier.f90:13
task_reconstruct
subroutine task_reconstruct(context)
Reconstructs the equilibrium.
Definition: v3fita.f:329
init_limiter_signals
subroutine init_limiter_signals(context, signals_created)
Constructs limiter signals.
Definition: v3fita.f:1648
file_opts::is_absolute_path
logical function is_absolute_path(path)
Queries if the a path is absoulte.
Definition: file_opts.f:309
integration_path
Module is part of the LIBSTELL. This modules contains code to define and integrate along an arbitray ...
Definition: integration_path.f:12
coordinate_utilities::cood_utils_test
logical function, public cood_utils_test()
Coordinate utilities unit test function.
Definition: coordinate_utilities.f:157
file_opts::delete_directory
subroutine delete_directory(directory_source, error)
Deletes a directory.
Definition: file_opts.f:279
thscte_dot
Module for opening and reading a 'thscte.' file. The file format for these files are documented in Th...
Definition: thscte_dot.f:73
init_vacuum_equilibrium
subroutine init_vacuum_equilibrium(context, force_solve)
Initialize a vacuum equilibrium.
Definition: v3fita.f:904
v3fit_context::v3fit_context_write
subroutine v3fit_context_write(this)
Write the v3fit_context_class out to disk.
Definition: v3fit_context.f:455
line_segment::line_seg_test
logical function, public line_seg_test()
Line segment unit test function.
Definition: line_segment.f:344
init_thomson_signals
subroutine init_thomson_signals(context, signals_created)
Constructs thomson signals.
Definition: v3fita.f:1405
v3fit_context::v3fit_context_resize
subroutine v3fit_context_resize(this)
Resize the arrays.
Definition: v3fit_context.f:277
data_parameters::process_dir
character(len=dir_prefix_len) function process_dir(rank)
Generates a directory name.
Definition: data_parameters.f:97
init_combination_signals
subroutine init_combination_signals(context, signals_created)
Constructs signal::combination_class signals.
Definition: v3fita.f:1891
v3fit_context::v3fit_context_destruct
subroutine v3fit_context_destruct(this)
Deconstruct a v3fit_context_class object.
Definition: v3fit_context.f:163
ipch_dot::ipch_dot_read
subroutine ipch_dot_read(ipch_file, signals, signals_created, observed, sigma, weight, first_index, use_polarimetry)
Read an interferometry/polarimetry diagnostic dot file.
Definition: ipch_dot.f:132
task_v3post
subroutine task_v3post(context)
Solves the equilibrium and calculates the modeled signals.
Definition: v3fita.f:227
file_opts
Contains cross platform routines for manipulating files on the file system. Defines a functions to mo...
Definition: file_opts.f:13
ece_dot::ece_dot_read
subroutine ece_dot_read(ece_file, signals, signals_created, observed, sigma, weight, first_index)
Read an ECE diagnostic dot file.
Definition: ece_dot.f:106
init_reconstruction
subroutine init_reconstruction(context)
Initialize Reconstruction.
Definition: v3fita.f:2103
init_vmec_equilibrium
subroutine init_vmec_equilibrium(context, force_solve)
Initialize a vmec equilibrium.
Definition: v3fita.f:822
sxrem_ratio_dot::sxrem_ratio_dot_read
subroutine sxrem_ratio_dot_read(sxrem_ratio_file, signals, signals_created, observed, sigma, weight, first_index)
Read an sxrem ratio dot file.
Definition: sxrem_ratio_dot.f:105
limiter_grid::limiter_grid_class
Base class representing a limiter signal.
Definition: limiter_grid.f:47
limiter_grid
Defines the base class of the type limiter_grid_class.
Definition: limiter_grid.f:13
init_prior_signals
subroutine init_prior_signals(context, signals_created)
Constructs prior signals.
Definition: v3fita.f:1762
integration_path::path_test
logical function path_test()
Path unit test function.
Definition: integration_path.f:851
file_opts::build_path
character(len=path_length) function build_path(path, name)
Builds a path.
Definition: file_opts.f:382
v3fit_context::v3fit_context_init_data
subroutine v3fit_context_init_data(this, eq_steps)
Initialize the dimensions and variables of the result file.
Definition: v3fit_context.f:808
profiler::profiler_construct
subroutine profiler_construct()
Construct a profiler.
Definition: profiler.f:67
init_signals
subroutine init_signals(context)
Constructs the signal array.
Definition: v3fita.f:1048
ipch_dot
Module for opening and reading a 'ipch.' file. The file format for these files are documented in Inte...
Definition: ipch_dot.f:100
file_opts::create_directory
subroutine create_directory(directory_source, error)
Makes a directory.
Definition: file_opts.f:249
data_parameters
This modules contains parameters used by equilibrium models.
Definition: data_parameters.f:10
line_segment
Module is part of the LIBSTELL. This module contains code to create a profile constructed of line sig...
Definition: line_segment.f:13
extcurz::extcurz_class
Base class representing a external Z currents enclosed by an inegration loop.
Definition: extcurz.f:35
task_equilibrium
subroutine task_equilibrium(context)
Solves the equilibrium.
Definition: v3fita.f:172
init_intpol_signals
subroutine init_intpol_signals(context, signals_created)
Constructs intpol signals.
Definition: v3fita.f:1340
init_sxrem_signals
subroutine init_sxrem_signals(context, signals_created)
Constructs sxrem signals.
Definition: v3fita.f:1285
functions
This module containes functions used by the profiles.
Definition: functions.f:10
init_parameters
subroutine init_parameters(context)
Constructs the v3fit_context::v3fit_context_class::params array.
Definition: v3fita.f:1975
profiler::profiler_write
subroutine profiler_write(iou)
Write out the profiled data to an output file.
Definition: profiler.f:336
v3fit_context::v3fit_context_class
Base class representing a v3fit context. This contains all memory needed to operate v3fit.
Definition: v3fit_context.f:30
v3fit_context::v3fit_context_construct
type(v3fit_context_class) function, pointer v3fit_context_construct(cl_parser)
Construct a v3fit_context_class object.
Definition: v3fit_context.f:130
v3fit_context::v3fit_context_write_step_data
subroutine v3fit_context_write_step_data(this, first_step, eq_steps)
Write step data to the defined variables.
Definition: v3fit_context.f:1123
prior_gaussian
Defines the base class of the type prior_gaussian_class. This class implements priors of the type.
Definition: prior_gaussian.f:17
limiter_iso_t::limiter_iso_class
Base class representing a limiter_iso function.
Definition: limiter_iso_T.f:43
init_gaussian_process
subroutine init_gaussian_process(context)
Constructs the v3fit_context::v3fit_context_class::gp array.
Definition: v3fita.f:2061
fourier::test_fourier
logical function test_fourier()
Test fourier routines.
Definition: fourier.f90:1050
cleanup_parallelism
subroutine cleanup_parallelism(context)
Cleanup Parallelism.
Definition: v3fita.f:2364
ece_dot
Module for opening and reading a 'ece.' file. The file format for these files are documented in ECE D...
Definition: ece_dot.f:75
task_child_equilibrium
subroutine task_child_equilibrium(context)
Child task to compute equilibrium.
Definition: v3fita.f:601
v3fit
program v3fit
V3FIT is a code for reconstructing 3D equilibria. V3FIT determined the optimal parameters of a 3D equ...
Definition: v3fita.f:47
limiter_iso_t
Defines the base class of the type limiter_iso.
Definition: limiter_iso_T.f:10