55 REAL (rprec) :: time_start, time_end
68 CALL second0(time_start)
75 IF (commandline_parser_is_flag_set(context%cl_parser,
79 CALL v3fit_input_read_namelist(
80 & commandline_parser_get_string(context%cl_parser,
89 CALL mpi_comm_rank(context%equilibrium_comm, eq_rank, error)
92 IF (eq_rank .eq. 0)
THEN
95 WRITE (context%recout_iou,*)
' my_task is ', my_task
96 WRITE (context%runlog_iou,*)
' my_task is ', my_task
104 CASE (
'v3post',
'vmec_v3post')
107 CASE (
'reconstruct',
'reconstruct_a1')
113 CASE (
'child_equilibrium')
120 CALL err_fatal(
'Unkown task type ' // trim(my_task))
125 IF (eq_rank .eq. 0)
THEN
128 CALL mpi_comm_rank(context%reconstruction_comm, recon_rank,
132 IF (recon_rank .eq. 0)
THEN
133 CALL second0(time_end)
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
153 CALL mpi_finalize(error)
156 1000
FORMAT(
' *** V3FIT run time (s): ',f10.2)
177 TYPE (v3fit_context_class),
INTENT(inout) :: context
183 REAL (rprec) :: start_time
186 start_time = profiler_get_start_time()
190 WRITE (*,*)
' *** In task equilibrium'
191 WRITE (context%runlog_iou,*)
' *** In task equilibrium'
197 WRITE (*,*)
' *** Starting initial equilibrium convergence (eq)'
198 WRITE (context%runlog_iou,*)
' *** Starting initial ' //
199 &
'equilibrium convergence'
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)
208 converged = model_converge(context%model, eq_steps,
209 & context%runlog_iou, 0,
'All')
214 CALL profiler_set_stop_time(
'task_equilibrium', start_time)
232 TYPE (v3fit_context_class),
INTENT(inout) :: context
239 INTEGER :: last_para_signal
241 REAL (rprec) :: start_time
244 REAL (rprec),
DIMENSION(4),
PARAMETER :: dummy_value = 0.0
247 start_time = profiler_get_start_time()
251 WRITE (*,*)
' *** In task v3post'
252 WRITE (context%runlog_iou,*)
' *** In task v3post'
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,
266 & context%equilibrium_comm,
'All')
267 CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0,
268 & context%equilibrium_comm, error)
274 DO i = 1,
SIZE(context%gp)
275 CALL gaussp_set_profile(context%gp(i)%p, context%model)
279 WRITE (*,*)
' *** Calculating modeled signals'
280 WRITE (context%runlog_iou,*)
' *** Calculating modeled signals'
287 IF (context%combination_index .eq. -1)
THEN
288 last_para_signal =
SIZE(context%signals)
290 last_para_signal = context%combination_index - 1
298 DO i = 1, last_para_signal
299 g2 = g2 + signal_get_g2(context%signals(i)%p,
300 & context%model, .false., dummy_value)
304 DO i = last_para_signal + 1,
SIZE(context%signals)
305 g2 = g2 + context%signals(i)%p%get_g2(context%model, .false.,
310 WRITE (context%runlog_iou,1000) g2
314 CALL profiler_set_stop_time(
'task_v3post', start_time)
316 1000
FORMAT(
'g^2 = ',es12.5)
334 TYPE (v3fit_context_class),
INTENT(inout) :: context
338 INTEGER :: restart_step
340 LOGICAL :: eq_converged
341 CHARACTER (len=dir_prefix_len) :: directory
342 LOGICAL :: write_input
343 REAL (rprec) :: start_time
346 start_time = profiler_get_start_time()
348 WRITE (*,*)
' *** In task reconstruct'
351 write_input = commandline_parser_is_flag_set(context%cl_parser,
360 IF (commandline_parser_is_flag_set(context%cl_parser,
362 WRITE (*,*)
' *** Restarting reconstruction '
368 CALL mpi_bcast(eq_steps, 1, mpi_integer, 0,
369 & context%reconstruction_comm, error)
373 WRITE (*,*)
' *** Starting initial equilibrium ' //
374 &
'convergence (rec)'
375 WRITE (context%runlog_iou,*)
' *** Starting initial ' //
376 &
'equilibrium convergence'
379 eq_converged = reconstruction_eval_e(context%recon,
381 & context%model, context%gp,
383 & context%runlog_iou,
385 & context%equilibrium_comm)
389 CALL assert(eq_converged,
390 &
'task_reconstruct no inital convergence')
391 CALL reconstruction_eval_f(context%recon,
392 & context%derived_params,
396 DO i = 1,
SIZE(context%signals)
397 context%recon%last_values(:,i) =
398 & context%signals(i)%p%modeled
404 CALL mpi_bcast(eq_steps, 1, mpi_integer, 0,
405 & context%reconstruction_comm, error)
410 WRITE (*,1003) reconstruction_get_g2(context%recon)
412 WRITE (context%runlog_iou,*)
413 WRITE (context%runlog_iou,1002) 0
414 WRITE (context%runlog_iou,1003)
415 & reconstruction_get_g2(context%recon)
417 IF (write_input)
THEN
418 CALL equilibrium_write_input(context%model%equilibrium, 0)
426 DO i = restart_step,
SIZE(context%recon%e, 2) - 1
430 WRITE (context%runlog_iou,*)
431 WRITE (context%runlog_iou,1002) i
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,
439 & context%reconstruction_comm,
440 & context%equilibrium_comm))
THEN
446 CALL reconstruction_eval_sem(context%recon,
449 & context%derived_params)
450 CALL reconstruction_write_step(context%recon,
451 & context%runlog_iou)
455 IF (write_input)
THEN
456 CALL equilibrium_write_input(context%model%equilibrium,
461 IF (reconstruction_get_dg2(context%recon) .le.
462 & context%recon_stop .or.
463 & reconstruction_get_g2(context%recon) .eq. 0.0)
THEN
465 WRITE (*,*)
' *** Equilibrium reconstructed'
466 WRITE (*,1000) context%recon_stop,
467 & reconstruction_get_dg2(context%recon)
469 WRITE (context%runlog_iou,*)
470 WRITE (context%runlog_iou,*)
' *** Equilibrium ' //
472 WRITE (context%runlog_iou,1000) context%recon_stop,
473 & reconstruction_get_dg2(context%recon)
477 IF (i .eq.
SIZE(context%recon%e, 2) - 1)
THEN
481 WRITE (context%runlog_iou,*)
482 WRITE (context%runlog_iou,1001) i
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)
496 CALL profiler_set_stop_time(
'task_reconstruct', start_time)
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)
516 TYPE (v3fit_context_class),
INTENT(inout) :: context
521 INTEGER :: signals_created
525 REAL (rprec) :: step_use
538 CALL mpi_bcast(eq_steps, 1, mpi_integer, 0,
539 & context%reconstruction_comm, error)
542 CALL mpi_bcast(task, 1, mpi_integer, 0,
543 & context%reconstruction_comm, error)
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)
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)
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)
579 CALL mpi_bcast(mpi_quit_task, 1, mpi_integer, 0,
580 & context%equilibrium_comm, error)
584 WRITE (*,*)
'Unknown MPI task.'
604 TYPE (v3fit_context_class),
INTENT(inout) :: context
617 CALL mpi_bcast(task, 1, mpi_integer, 0,
618 & context%equilibrium_comm, error)
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,
627 & context%runlog_iou,
628 & context%equilibrium_comm,
629 & context%model%state_flags)
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,
636 & context%equilibrium_comm)
642 WRITE (*,*)
'Unknown MPI task.'
664 TYPE (v3fit_context_class),
INTENT(inout) :: context
667 REAL (rprec) :: start_time
670 start_time = profiler_get_start_time()
672 WRITE (*,*)
' *** In task gsq_on_grid'
698 LOGICAL :: tests_failed
701 WRITE (*,*)
' *** In task unit_test'
702 tests_failed = .false.
706 WRITE (*,*)
'Line Segment Test Passed'
708 tests_failed = .true.
713 WRITE (*,*)
'Coord Utilites Test Passed'
715 tests_failed = .true.
719 IF (function_test())
THEN
720 WRITE (*,*)
'Functions Test Passed'
722 tests_failed = .true.
727 WRITE (*,*)
'Integration Path Test Passed'
729 tests_failed = .true.
734 WRITE (*,*)
'SIESTA Fourier Test Passed'
736 tests_failed = .true.
741 IF (tests_failed)
THEN
765 TYPE (v3fit_context_class),
INTENT(inout) :: context
768 LOGICAL :: force_solve
771 REAL (rprec) :: start_time
774 start_time = profiler_get_start_time()
778 CALL mpi_comm_rank(context%equilibrium_comm, eq_rank, error)
781 WRITE (*,*)
' *** Constructing equilibrium model'
782 IF (eq_rank .eq. 0)
THEN
783 WRITE (context%runlog_iou,*)
' *** Constructing ' //
784 &
'equilibrium model'
788 & commandline_parser_is_flag_set(context%cl_parser,
'-force')
790 SELECT CASE (model_eq_type)
802 WRITE (*,*)
'Equilibrium not initialized. Please ' //
803 &
'specify a vaild equilibrium type.'
808 CALL profiler_set_stop_time(
'init_equilibrium', start_time)
828 TYPE (v3fit_context_class),
INTENT(inout) :: context
829 LOGICAL,
INTENT(in) :: force_solve
832 TYPE (pprofile_pointer),
DIMENSION(:),
POINTER :: sxr => null()
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
841 start_time = profiler_get_start_time()
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)),
847 & pp_sxrem_as_a(i,:),
848 & pp_sxrem_af_a(i,:))
852 IF (trim(emission_file) .ne.
'')
THEN
853 emission_func => emission_construct(emission_file)
856 state_flags = model_state_all_off
857 vmec => vmec_construct(vmec_nli_filename, vmec_wout_input,
858 & pprofile_construct(trim(pp_ne_ptype),
861 & pprofile_construct(trim(pp_te_ptype),
864 & pprofile_construct(trim(pp_ti_ptype),
867 & pprofile_construct(trim(pp_ze_ptype),
870 & phi_offset, z_offset, pol_rad_ratio,
871 & context%runlog_iou,
873 & context%equilibrium_comm,
874 & context%reconstruction_comm, state_flags)
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,
886 & sfactor_spec_fac(:minloc(sfactor_spec_imin, 1) - 1),
887 & soffset_spec_fac(:minloc(soffset_spec_imin, 1) - 1),
890 CALL profiler_set_stop_time(
'init_vmec_equilibrium', start_time)
910 TYPE (v3fit_context_class),
INTENT(inout) :: context
911 LOGICAL,
INTENT(in) :: force_solve
914 REAL (rprec) :: start_time
915 TYPE (emission_class),
POINTER :: emission_func => null()
918 start_time = profiler_get_start_time()
921 IF (trim(emission_file) .ne.
'')
THEN
922 emission_func => emission_construct(emission_file)
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,
930 & equilibrium_construct(vacuum_construct(vacuum_nli_filename,
931 & context%runlog_iou),
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),
939 CALL profiler_set_stop_time(
'init_vacuum_equilibrium', start_time)
959 TYPE (v3fit_context_class),
INTENT(inout) :: context
960 LOGICAL,
INTENT(in) :: force_solve
963 REAL (rprec) :: start_time
964 TYPE (pprofile_pointer),
DIMENSION(:),
POINTER :: sxr => null()
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
973 start_time = profiler_get_start_time()
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)),
979 & pp_sxrem_as_a(i,:),
980 & pp_sxrem_af_a(i,:))
984 IF (trim(emission_file) .ne.
'')
THEN
985 emission_func => emission_construct(emission_file)
988 state_flags = model_state_all_off
989 vmec => vmec_construct(vmec_nli_filename, vmec_wout_input,
990 & context%runlog_iou,
992 & context%equilibrium_comm,
993 & context%reconstruction_comm,
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,
1004 & pprofile_construct(trim(pp_te_ptype),
1005 & pp_te_b, pp_te_as,
1007 & pprofile_construct(trim(pp_ti_ptype),
1008 & pp_ti_b, pp_ti_as,
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,
1018 & vmec, state_flags)
1022 not_converged = trim(siesta_restart_filename) .eq.
''
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,
1030 & sxrem_te_a, sxrem_ratio_a, ece_resonance_range, coosig_wgts,
1032 & sfactor_spec_fac(:minloc(sfactor_spec_imin, 1) - 1),
1033 & soffset_spec_fac(:minloc(soffset_spec_imin, 1) - 1),
1036 CALL profiler_set_stop_time(
'init_siesta_equilibrium', start_time)
1053 TYPE (v3fit_context_class),
INTENT(inout) :: context
1056 INTEGER :: signals_created = 0
1059 REAL (rprec) :: start_time
1062 start_time = profiler_get_start_time()
1064 IF (na_s_desc .eq. 0)
THEN
1065 CALL err_warn(
'init_signals: no signals will be created. ' //
1070 WRITE (*,*)
' *** Constructing signals'
1071 WRITE (context%runlog_iou,*)
' *** Constructing signals'
1073 ALLOCATE(context%signals(na_s_desc))
1076 IF (mdsig_list_filename .ne.
'')
THEN
1081 IF (sxrch_dot_filename .ne.
'')
THEN
1086 IF (ipch_dot_filename .ne.
'')
THEN
1091 IF (thscte_dot_filename .ne.
'')
THEN
1096 IF (extcurz_s0 .ge. 0.0)
THEN
1101 IF (mse_dot_filename .ne.
'')
THEN
1106 IF (ece_dot_filename .ne.
'')
THEN
1114 IF (sxrem_ratio_dot_filename .ne.
'')
THEN
1122 CALL profiler_set_stop_time(
'init_signals', start_time)
1147 TYPE (v3fit_context_class),
INTENT(inout) :: context
1148 INTEGER,
INTENT(inout) :: signals_created
1151 INTEGER :: plasma_index
1152 LOGICAL :: use_point
1153 INTEGER :: mdsig_list_iou
1154 INTEGER :: mdsig_iou
1156 INTEGER :: mdsig_index
1157 CHARACTER (len=path_length) :: mdsig_path
1158 CHARACTER (len=path_length) :: mdsig_filename
1161 REAL (rprec) :: start_time
1164 start_time = profiler_get_start_time()
1170 #if defined(MPI_OPT)
1171 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, status)
1173 IF (mpi_rank .gt. 0)
THEN
1176 mdsig_list_filename =
build_path(
'..', mdsig_list_filename)
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)
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))
1200 & commandline_parser_get_real(context%cl_parser,
'-c',
1205 READ(mdsig_list_iou, *, iostat=status)
1206 & mdsig_index, mdsig_filename
1210 IF (status .ne. 0)
EXIT
1213 IF (signals_created + 1 .gt. v3fit_max_diagnostics)
THEN
1214 CALL err_fatal(
'init_magnetic_signals: created signals' //
1215 &
' exceeds v3fit_max_diagnostics')
1219 CALL cdf_open(mdsig_iou, trim(
build_path(mdsig_path,
1222 CALL assert_eq(0, status,
'init_magnetic_signals: ' //
1223 &
'failed to open ', trim(mdsig_path) //
1224 & trim(mdsig_filename))
1227 signals_created = signals_created + 1
1230 & mag_force_coil, mag_3d_a(signals_created),
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))
1240 CALL cdf_close(mdsig_iou)
1242 context%signals(signals_created)%p => magnetic_object
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)
1253 IF (.not.use_point .and.
1254 & magnetic_response_is_point(magnetic_object%response))
THEN
1256 CALL equilibrium_set_magnetic_cache(
1257 & context%model%equilibrium, any(mag_3d_a),
1258 & context%model%state_flags)
1263 IF (context%magnetic_index .eq. -1)
THEN
1264 context%magnetic_index = signals_created
1268 CLOSE(mdsig_list_iou)
1270 CALL profiler_set_stop_time(
'init_magnetic_signals', start_time)
1292 TYPE (v3fit_context_class),
INTENT(inout) :: context
1293 INTEGER,
INTENT(inout) :: signals_created
1298 REAL (rprec) :: start_time
1301 start_time = profiler_get_start_time()
1304 #if defined(MPI_OPT)
1305 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1307 IF (mpi_rank .gt. 0)
THEN
1310 sxrch_dot_filename =
build_path(
'..', sxrch_dot_filename)
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)
1322 & signals_created, sdo_data_a, sdo_sigma_a,
1323 & sdo_weight_a, context%sxrem_index)
1325 CALL profiler_set_stop_time(
'init_sxrem_signals', start_time)
1347 TYPE (v3fit_context_class),
INTENT(inout) :: context
1348 INTEGER,
INTENT(inout) :: signals_created
1351 LOGICAL :: use_polarimetry
1354 REAL (rprec) :: start_time
1357 start_time = profiler_get_start_time()
1360 #if defined(MPI_OPT)
1361 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1363 IF (mpi_rank .gt. 0)
THEN
1366 ipch_dot_filename =
build_path(
'..', ipch_dot_filename)
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)
1378 & signals_created, sdo_data_a, sdo_sigma_a,
1379 & sdo_weight_a, context%intpol_index,
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,
1387 & context%model%state_flags)
1390 CALL profiler_set_stop_time(
'init_intpol_signals', start_time)
1412 TYPE (v3fit_context_class),
INTENT(inout) :: context
1413 INTEGER,
INTENT(inout) :: signals_created
1418 REAL (rprec) :: start_time
1421 start_time = profiler_get_start_time()
1424 #if defined(MPI_OPT)
1425 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1427 IF (mpi_rank .gt. 0)
THEN
1430 thscte_dot_filename =
build_path(
'..', thscte_dot_filename)
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)
1443 & signals_created, sdo_data_a, sdo_sigma_a,
1444 & sdo_weight_a, context%thomson_index)
1446 CALL profiler_set_stop_time(
'init_thomson_signals', start_time)
1467 TYPE (v3fit_context_class),
INTENT(inout) :: context
1468 INTEGER,
INTENT(inout) :: signals_created
1473 class(signal_class),
POINTER :: extcurz_object
1474 REAL (rprec) :: start_time
1477 start_time = profiler_get_start_time()
1479 WRITE (*,*)
' *** Constructing extcurz signals'
1480 WRITE (context%runlog_iou,*)
' *** Constructing extcurz 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')
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))
1498 context%signals(signals_created + 1)%p => extcurz_object
1500 signals_created = signals_created + 1
1504 IF (context%extcurz_index .eq. -1)
THEN
1505 context%extcurz_index = signals_created
1508 CALL profiler_set_stop_time(
'init_extcurz_signals', start_time)
1530 TYPE (v3fit_context_class),
INTENT(inout) :: context
1531 INTEGER,
INTENT(inout) :: signals_created
1536 REAL (rprec) :: start_time
1539 start_time = profiler_get_start_time()
1542 #if defined(MPI_OPT)
1543 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1545 IF (mpi_rank .gt. 0)
THEN
1548 mse_dot_filename =
build_path(
'..', mse_dot_filename)
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)
1560 & signals_created, sdo_data_a, sdo_sigma_a,
1561 & sdo_weight_a, context%mse_index)
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,
1568 & context%model%state_flags)
1571 CALL profiler_set_stop_time(
'init_mse_signals', start_time)
1593 TYPE (v3fit_context_class),
INTENT(inout) :: context
1594 INTEGER,
INTENT(inout) :: signals_created
1599 REAL (rprec) :: start_time
1602 start_time = profiler_get_start_time()
1605 #if defined(MPI_OPT)
1606 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1608 IF (mpi_rank .gt. 0)
THEN
1611 ece_dot_filename =
build_path(
'..', ece_dot_filename)
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)
1623 & signals_created, sdo_data_a, sdo_sigma_a,
1624 & sdo_weight_a, context%ece_index)
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,
1631 & context%model%state_flags)
1634 CALL profiler_set_stop_time(
'init_ece_signals', start_time)
1649 USE stel_constants,
only: degree
1656 TYPE (v3fit_context_class),
INTENT(inout) :: context
1657 INTEGER,
INTENT(inout) :: signals_created
1661 CHARACTER (len=data_short_name_length) :: s_name =
''
1662 CHARACTER (len=data_name_length) :: l_name =
''
1664 class(signal_class),
POINTER :: limiter_obj
1665 REAL (rprec) :: start_time
1668 start_time = profiler_get_start_time()
1670 WRITE (*,*)
' *** Constructing limiter signals'
1671 WRITE (context%runlog_iou,*)
' *** Constructing limiter 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')
1682 IF (i .gt. v3fit_max_limiters)
THEN
1683 CALL err_fatal(
'init_limiter_signals: created signals' //
1684 &
' exceeds v3fit_max_limits')
1687 s_name = trim(signal_make_short_name(
'edge_lim_', i))
1688 WRITE (l_name,
'(a,i4)')
'edge limiter iso function ', i
1690 num_phi = max(min(v3fit_max_lif_size, n_phi_lif(i)), 1)
1693 IF (sdo_sigma_a(signals_created + 1) .eq. 0.0)
THEN
1694 sdo_sigma_a(signals_created + 1) = lif_sigma(i)
1699 & lif_phi_degree(i,1:num_phi)*degree,
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))
1709 context%signals(signals_created + 1)%p => limiter_obj
1711 signals_created = signals_created + 1
1715 IF (context%limiter_index .eq. -1)
THEN
1716 context%limiter_index = signals_created
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
1727 & lif_on_edge(n_lif + 1))
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))
1736 context%signals(signals_created + 1)%p => limiter_obj
1738 signals_created = signals_created + 1
1742 IF (context%limiter_index .eq. -1)
THEN
1743 context%limiter_index = signals_created
1748 CALL profiler_set_stop_time(
'init_limiter_signals', start_time)
1768 TYPE (v3fit_context_class),
INTENT(inout) :: context
1769 INTEGER,
INTENT(inout) :: signals_created
1772 class(signal_class),
POINTER :: prior_obj
1774 REAL (rprec) :: start_time
1777 start_time = profiler_get_start_time()
1779 WRITE (*,*)
' *** Constructing prior signals'
1780 WRITE (context%runlog_iou,*)
' *** Constructing prior 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')
1791 IF (i .gt. v3fit_max_priors)
THEN
1792 CALL err_fatal(
'init_prior_signals: created signals' //
1793 &
' exceeds v3fit_max_priors')
1797 & prior_param_name(i),
1798 & prior_indices(i,:))
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))
1809 context%signals(signals_created + 1)%p => prior_obj
1811 signals_created = signals_created + 1
1815 IF (context%prior_gaussian_index .eq. -1)
THEN
1816 context%prior_gaussian_index = signals_created
1821 CALL profiler_set_stop_time(
'init_prior_signals', start_time)
1843 TYPE (v3fit_context_class),
INTENT(inout) :: context
1844 INTEGER,
INTENT(inout) :: signals_created
1849 REAL (rprec) :: start_time
1852 start_time = profiler_get_start_time()
1855 #if defined(MPI_OPT)
1856 CALL mpi_comm_rank(context%reconstruction_comm, mpi_rank, error)
1858 IF (mpi_rank .gt. 0)
THEN
1861 sxrem_ratio_dot_filename =
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)
1874 & context%signals, signals_created,
1875 & sdo_data_a, sdo_sigma_a, sdo_weight_a,
1876 & context%sxrem_ratio_index)
1878 CALL profiler_set_stop_time(
'init_sxrem_ratio_signals',
1896 TYPE (v3fit_context_class),
INTENT(inout) :: context
1897 INTEGER,
INTENT(inout) :: signals_created
1902 class(combination_class),
POINTER :: combination_signal
1903 REAL (rprec) :: start_time
1906 start_time = profiler_get_start_time()
1908 WRITE (*,*)
' *** Constructing combination signals'
1909 WRITE (context%runlog_iou,*)
' *** Constructing ' //
1910 &
'combination 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')
1920 IF (i .gt. v3fit_max_combinations)
THEN
1921 CALL err_fatal(
'init_combination_signals: created ' //
1922 &
'signals exceeds v3fit_max_combinations')
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')
1936 CALL combination_set_signal(combination_signal,
1937 & context%signals(coosig_indices(i,j))%p,
1938 & coosig_coeff(i,j), j)
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))
1949 context%signals(signals_created + 1)%p => combination_signal
1951 signals_created = signals_created + 1
1955 IF (context%combination_index .eq. -1)
THEN
1956 context%combination_index = signals_created
1961 CALL profiler_set_stop_time(
'init_combination_signals',
1980 TYPE (v3fit_context_class),
INTENT(inout) :: context
1984 INTEGER,
DIMENSION(:,:),
ALLOCATABLE :: indices
1985 REAL (rprec) :: start_time
1988 start_time = profiler_get_start_time()
1990 WRITE (*,*)
' *** Constructing reconstruction parameters'
1991 WRITE (context%runlog_iou,*)
' *** Constructing ' //
1992 &
'reconstruction parameters'
1995 ALLOCATE(context%params(n_rp))
1997 context%params(i)%p => param_construct(context%model,
2002 & rp_range_type(i,:),
2003 & rp_range_index(i,:,:),
2004 & rp_range_value(i,:),
2005 &
SIZE(context%signals),
2010 ALLOCATE(context%derived_params(n_dp))
2012 context%derived_params(i)%p => param_construct(context%model,
2020 ALLOCATE(context%locks(n_lp))
2024 DO j = 1,
SIZE(lp_sets,2)
2025 IF (lp_sets(i,j) .eq.
'')
THEN
2029 j = min(j - 1,
SIZE(lp_sets,2))
2031 CALL assert(j .ne. 0,
'Locking parameters require at ' //
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)
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)
2048 CALL profiler_set_stop_time(
'init_parameters', start_time)
2066 TYPE (v3fit_context_class),
INTENT(inout) :: context
2071 REAL (rprec) :: start_time
2074 start_time = profiler_get_start_time()
2076 ALLOCATE(context%gp(n_gp))
2079 context%gp(i)%p => gaussp_construct(context%model,
2082 & gp_model_index(i),
2083 & gp_param_vrnc(i,:),
2085 & gp_cholesky_fact(i))
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)
2093 CALL profiler_set_stop_time(
'init_gaussian_process', start_time)
2108 TYPE (v3fit_context_class),
INTENT(inout) :: context
2111 INTEGER :: last_para_signal
2112 REAL (rprec) :: start_time
2115 start_time = profiler_get_start_time()
2120 IF (context%combination_index .eq. -1)
THEN
2121 last_para_signal =
SIZE(context%signals)
2123 last_para_signal = context%combination_index - 1
2127 use_central_diff = use_central_diff .or.
2128 & commandline_parser_is_flag_set(context%cl_parser,
'-c_diff')
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,
2140 CALL profiler_set_stop_time(
'init_reconstruction', start_time)
2175 #if defined(MPI_OPT)
2182 TYPE (v3fit_context_class),
INTENT(inout) :: context
2185 INTEGER :: num_threads
2186 #if defined(MPI_OPT)
2188 INTEGER :: num_recon_processes
2189 INTEGER :: num_eq_processes
2195 INTEGER :: recon_rank
2196 CHARACTER (len=dir_prefix_len) :: directory_name
2197 CHARACTER (len=path_length) :: temp_string
2199 REAL (rprec) :: start_time
2202 start_time = profiler_get_start_time()
2204 WRITE (*,*)
' *** Configuring parallelism.'
2222 #if defined (MPI_OPT)
2224 CALL mpi_comm_size(mpi_comm_world, mpi_size, error)
2225 CALL mpi_comm_rank(mpi_comm_world, mpi_rank, error)
2230 context%global_comm = mpi_comm_world
2234 SELECT CASE (my_task)
2236 CASE (
'reconstruct',
'reconstruct_a1')
2237 IF (commandline_parser_is_flag_set(context%cl_parser,
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
2245 num_recon_processes = n_rp
2246 num_eq_processes = mpi_size/n_rp
2247 & + mod(mpi_size, n_rp)
2250 CASE (
'v3post',
'vmec_v3post',
'equilibrium')
2251 num_recon_processes = 1
2252 num_eq_processes = mpi_size
2255 num_recon_processes = 1
2256 IF (mpi_rank .gt. 0)
THEN
2257 stop
'MPI is only supported for the reconstruct task.'
2266 temp_string = commandline_parser_get_string(context%cl_parser,
2269 IF (mpi_rank .eq. 0)
THEN
2270 DO i = 2, num_recon_processes
2284 & vmec_nli_filename), error)
2289 & vmec_wout_input), error)
2294 & vacuum_nli_filename), error)
2299 & siesta_nli_filename), error)
2304 & siesta_restart_filename),
2311 CALL mpi_barrier(context%global_comm, error)
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)
2321 IF (eq_rank .gt. 0)
THEN
2322 CALL mpi_comm_split(context%global_comm, mpi_undefined, 0,
2323 & context%reconstruction_comm, error)
2325 CALL mpi_comm_split(context%global_comm, 0, recon_rank,
2326 & context%reconstruction_comm, error)
2330 IF (recon_rank .gt. 0)
THEN
2334 IF (mpi_rank .gt. 0)
THEN
2335 IF (eq_rank .gt. 0)
THEN
2336 my_task =
'child_equilibrium'
2338 my_task =
'child_recon'
2342 WRITE (*,1001) mpi_rank + 1, recon_rank + 1, eq_rank + 1,
2345 CALL mpi_barrier(context%global_comm, error)
2348 CALL profiler_set_stop_time(
'config_parallelism', start_time)
2350 1000
FORMAT(
'Using ',i4,
' threads.')
2351 1001
FORMAT(
'Configured global process ',i5,
', recon process ',i5,
2352 &
', eq process ',i5,
' task ',a)
2366 #if defined(MPI_OPT)
2373 TYPE (v3fit_context_class),
INTENT(inout) :: context
2376 #if defined(MPI_OPT)
2380 INTEGER :: recon_size
2384 IF (context%reconstruction_comm .eq. mpi_comm_null)
THEN
2388 CALL mpi_comm_rank(context%global_comm, mpi_rank, error)
2390 CALL mpi_comm_size(context%reconstruction_comm, recon_size,
2392 CALL mpi_barrier(context%reconstruction_comm, error)
2395 IF (mpi_rank .eq. 0)
THEN
2398 DO i = 2, recon_size