Merge pull request #206 from PrincetonUniversity/bugfix/BEAMS3D_minor

Bugfix/beams3 d minor
PrincetonUniversity · Oct 26, 2023 · af17e3a · af17e3a
2 parents d374529 + 631f90e
commit af17e3a
Show file tree

Hide file tree

Showing 7 changed files with 317 additions and 255 deletions.
diff --git a/BEAMS3D/BEAMS3D.dep b/BEAMS3D/BEAMS3D.dep
@@ -366,6 +366,7 @@ fidasim_input_mod.o: \
       ../../LIBSTELL/$(LOCTYPE)/wall_mod.o \
       ../../LIBSTELL/$(LOCTYPE)/mpi_params.o \
       ../../LIBSTELL/$(LOCTYPE)/mpi_inc.o \
+      beams3d_physics_mod.o \
       beams3d_runtime.o \
       beams3d_lines.o \
       beams3d_write_parhdf5.o \

diff --git a/BEAMS3D/Sources/beams3d_diagnostics.f90 b/BEAMS3D/Sources/beams3d_diagnostics.f90
@@ -107,8 +107,8 @@ SUBROUTINE beams3d_diagnostics
          CALL MPI_REDUCE(shine_port,    shine_port,    nbeams,                  MPI_DOUBLE_PRECISION, MPI_SUM, master, MPI_COMM_BEAMS, ierr_mpi)
          CALL MPI_REDUCE(norbit,        norbit,        nbeams,                  MPI_REAL,          MPI_SUM, master, MPI_COMM_BEAMS, ierr_mpi)
          CALL MPI_REDUCE(nlost,         nlost,         nbeams,                  MPI_REAL,          MPI_SUM, master, MPI_COMM_BEAMS, ierr_mpi)
-         CALL MPI_REDUCE(nlost,         tlow,          nbeams,                  MPI_REAL,          MPI_MIN, master, MPI_COMM_BEAMS, ierr_mpi)
-         CALL MPI_REDUCE(nlost,         thigh,         nbeams,                  MPI_REAL,          MPI_MAX, master, MPI_COMM_BEAMS, ierr_mpi)
+         CALL MPI_REDUCE(tlow,          tlow,          nbeams,                  MPI_REAL,          MPI_MIN, master, MPI_COMM_BEAMS, ierr_mpi)
+         CALL MPI_REDUCE(thigh,         thigh,         nbeams,                  MPI_REAL,          MPI_MAX, master, MPI_COMM_BEAMS, ierr_mpi)
       END IF
 #endif
 

diff --git a/BEAMS3D/Sources/beams3d_follow.f90 b/BEAMS3D/Sources/beams3d_follow.f90
@@ -302,11 +302,14 @@ SUBROUTINE beams3d_follow
     DEALLOCATE(itemp)
     IF (ALLOCATED(mnum)) DEALLOCATE(mnum)
     IF (ALLOCATED(moffsets)) DEALLOCATE(moffsets)
-    IF (ALLOCATED(t_last)) DEALLOCATE(t_last)
     CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
     IF (ierr_mpi /= 0) CALL handle_err(MPI_BARRIER_ERR, 'beams3d_follow', ierr_mpi)
 #endif
 
+    ! Adjust T_END back to values of T_last
+    t_end(mystart:myend) = t_last(mystart:myend)
+    IF (ALLOCATED(t_last)) DEALLOCATE(t_last)
+
     RETURN
     !-----------------------------------------------------------------------
     !     End Subroutine

diff --git a/BEAMS3D/Sources/beams3d_init.f90 b/BEAMS3D/Sources/beams3d_init.f90
@@ -423,7 +423,7 @@ SUBROUTINE beams3d_init
       ! Load vessel if not done already vessel
       IF (lvessel .and. (.not. lwall_loaded)) THEN
       IF (lverb) THEN
-         WRITE(6,'(A)') 'Loading Vessel!'
+         WRITE(6,'(A)') '----- Loading wall data -----'
       END IF
          CALL wall_load_txt(TRIM(vessel_string),ier,.false.,MPI_COMM_BEAMS)
          IF (lverb) THEN

diff --git a/BEAMS3D/Sources/beams3d_interface_mod.f90 b/BEAMS3D/Sources/beams3d_interface_mod.f90
@@ -78,12 +78,15 @@ SUBROUTINE beams3d_cleanup
       CALL beams3d_free(MPI_COMM_SHARMEM)
       IF (lvessel) CALL wall_free(ier,MPI_COMM_BEAMS)
 #if defined(MPI_OPT)
-      CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
+      ierr_mpi = 0; CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
       IF (ierr_mpi /= 0) CALL handle_err(MPI_BARRIER_ERR, 'beams3d_main', ierr_mpi)
-      ierr_mpi=0
-      CALL MPI_INFO_FREE(mpi_info_beams3d, ierr_mpi)
-      CALL MPI_FINALIZE(ierr_mpi)
-      IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_FINE_ERR, 'beams3d_main', ierr_mpi)
+      ierr_mpi = 0; CALL MPI_INFO_FREE(mpi_info_beams3d, ierr_mpi)
+      ierr_mpi = 0; CALL MPI_BARRIER(MPI_COMM_BEAMS, ierr_mpi)
+      IF (ierr_mpi /= 0) CALL handle_err(MPI_BARRIER_ERR, 'beams3d_main', ierr_mpi)
+      ierr_mpi = 0; CALL MPI_COMM_FREE(MPI_COMM_SHARMEM, ierr_mpi)
+      ierr_mpi = 0; CALL MPI_COMM_FREE(MPI_COMM_BEAMS, ierr_mpi)
+      ierr_mpi = 0; CALL MPI_FINALIZE(ierr_mpi)
+      !IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_FINE_ERR, 'beams3d_main', ierr_mpi)
 #endif
       IF (lverb) WRITE(6, '(A)') '----- BEAMS3D DONE -----'
       RETURN

diff --git a/BEAMS3D/Sources/beams3d_physics_mod.f90 b/BEAMS3D/Sources/beams3d_physics_mod.f90
@@ -1811,6 +1811,78 @@ SUBROUTINE beams3d_SFLX(q,S)
 
       END SUBROUTINE beams3d_SFLX
 
+      !-----------------------------------------------------------------
+      !     Function:      beams3d_SUFLX
+      !     Authors:       S. Lazerson ([email protected])
+      !     Date:          09/28/2023
+      !     Description:   Returns normalized toroidal flux and
+      !                    poloidal angle at a given point in space.
+      !-----------------------------------------------------------------
+      SUBROUTINE beams3d_SUFLX(q,S,U)
+         !--------------------------------------------------------------
+         !     Input Parameters
+         !          q    (q(1),q(2),q(3)) = (R,phi,Z)
+         !          S    Backbround grid normalized flux
+         !          U    Backbround grid poloidal angle
+         !--------------------------------------------------------------
+         IMPLICIT NONE
+         DOUBLE PRECISION, INTENT(inout) :: q(3)
+         DOUBLE PRECISION, INTENT(out) :: S
+         DOUBLE PRECISION, INTENT(out) :: U
+
+         !--------------------------------------------------------------
+         !     Local Variables
+         !        r_temp     Helpers (r,phi,z)
+         !        i,j,k      Spline Grid indicies
+         !        xparam     Spline subgrid factor [0,1] (yparam,zparam)
+         !        ict        Spline output control
+         !        fval       Spline output array
+         !--------------------------------------------------------------
+         DOUBLE PRECISION :: r_temp, z_temp, phi_temp
+         ! For splines
+         INTEGER :: i,j,k
+         REAL*8 :: xparam, yparam, zparam
+         INTEGER, parameter :: ict(8)=(/1,0,0,0,0,0,0,0/)
+         REAL*8 :: fval(1)
+
+         !--------------------------------------------------------------
+         !     Begin Subroutine
+         !--------------------------------------------------------------
+
+         ! Setup position in a vll arrays
+         r_temp   = q(1)
+         phi_temp = MODULO(q(2), phimax)
+         IF (phi_temp < 0) phi_temp = phi_temp + phimax
+         z_temp   = q(3)
+
+         ! Initialize values
+         S = 2; U = 0
+
+         ! Check that we're inside the domain then proceed
+         IF ((r_temp >= rmin-eps1) .and. (r_temp <= rmax+eps1) .and. &
+             (phi_temp >= phimin-eps2) .and. (phi_temp <= phimax+eps2) .and. &
+             (z_temp >= zmin-eps3) .and. (z_temp <= zmax+eps3)) THEN
+            i = MIN(MAX(COUNT(raxis < r_temp),1),nr-1)
+            j = MIN(MAX(COUNT(phiaxis < phi_temp),1),nphi-1)
+            k = MIN(MAX(COUNT(zaxis < z_temp),1),nz-1)
+            xparam = (r_temp - raxis(i)) * hri(i)
+            yparam = (phi_temp - phiaxis(j)) * hpi(j)
+            zparam = (z_temp - zaxis(k)) * hzi(k)
+            ! Evaluate the Splines
+            CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
+                            hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
+                            S4D(1,1,1,1),nr,nphi,nz)
+            S = max(fval(1),zero)
+            CALL R8HERM3FCN(ict,1,1,fval,i,j,k,xparam,yparam,zparam,&
+                            hr(i),hri(i),hp(j),hpi(j),hz(k),hzi(k),&
+                            U4D(1,1,1,1),nr,nphi,nz)
+            U = fval(1)
+         END IF
+
+         RETURN
+
+      END SUBROUTINE beams3d_SUFLX
+
       !-----------------------------------------------------------------
       !     Function:      beams3d_suv2rzp
       !     Authors:       S. Lazerson ([email protected])