diff --git a/Utilities/pythontools/py_spec/ci/test.py b/Utilities/pythontools/py_spec/ci/test.py index 57708c5d..94b1411e 100755 --- a/Utilities/pythontools/py_spec/ci/test.py +++ b/Utilities/pythontools/py_spec/ci/test.py @@ -23,6 +23,8 @@ def compare(data, reference, localtol=1e-6, action='ERR'): """ global match for key, value in vars(data).items(): + if key in ['Lconstraint', 'LBeltrami']: + continue if isinstance(value, SPECout): # recurse data (csmiet: I'm not the biggest fan of this recursion...) print('------------------') print('Elements in '+key) diff --git a/ci/G1V03L2Fi/G1V03L2Fi.001.sp b/ci/G1V03L2Fi/G1V03L2Fi.001.sp index 31e423fc..47024fc5 100644 --- a/ci/G1V03L2Fi/G1V03L2Fi.001.sp +++ b/ci/G1V03L2Fi/G1V03L2Fi.001.sp @@ -63,7 +63,7 @@ Vns(0,2) = 0.000000000000000E+00 Bns(0,2) = 0.000000000000000E+00 Vnc(0, Mregular = -1 / &locallist - LBeltrami = 2 + LBeltrami = 4 Linitgues = 1 / &globallist diff --git a/src/dforce.f90 b/src/dforce.f90 index 05f59b2c..ccc7da32 100644 --- a/src/dforce.f90 +++ b/src/dforce.f90 @@ -98,14 +98,14 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) use numerical, only : logtolerance - use fileunits, only : ounit + use fileunits, only : ounit, munit use inputlist, only : Wmacros, Wdforce, Nvol, Ntor, Lrad, Igeometry, & epsilon, & Lconstraint, Lcheck, dRZ, & Lextrap, & mupftol, & - Lfreebound, LHmatrix + Lfreebound, LHmatrix, gamma, pscale, adiabatic use cputiming, only : Tdforce @@ -134,7 +134,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) LocalConstraint, xoffset, & solution, IPdtdPf, & IsMyVolume, IsMyVolumeValue, WhichCpuID, & - ext ! For outputing Lcheck = 6 test + ext, & ! For outputing Lcheck = 6 test + BetaTotal, vvolume !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -164,6 +165,9 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) LOGICAL :: LComputeAxis, dfp100_logical + REAL :: press, voltotal + REAL :: betavol(1:Mvol) + #ifdef DEBUG INTEGER :: isymdiff REAL :: dvol(-1:+1), evolume, imupf(1:2,-2:2), lfactor @@ -239,8 +243,13 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) ! Mvol-1 surface current plus 1 poloidal linking current constraints Ndofgl = Mvol else - ! Mvol-1 surface current constraints - Ndofgl = Mvol-1 + ! add an additional constraint to make the total pflux = 0 + if(Igeometry.eq.1) then + Ndofgl = Mvol + else + ! Mvol-1 surface current constraints + Ndofgl = Mvol-1 + endif endif SALLOCATE( Fvec, (1:Ndofgl), zero ) @@ -257,6 +266,11 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) ! one step Newton's method dpflux(2:Mvol) = dpflux(2:Mvol) - dpfluxout(1:Mvol-1) + + if(Igeometry.eq.1) then + dpflux(1) = dpflux(1) - dpfluxout(Mvol) + endif + if( Lfreebound.eq.1 ) then dtflux(Mvol) = dtflux(Mvol ) - dpfluxout(Mvol ) endif @@ -305,6 +319,38 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) enddo ! end of do vvol = 1, Mvol + !add an additional constraint to make the total pflux 0 + if(Igeometry.eq.1) then + + vvol = 1 + + WCALL(dforce, IsMyVolume, (vvol)) + + if( IsMyVolumeValue .EQ. 0 ) then + + else if( IsMyVolumeValue .EQ. -1) then + FATAL(dforce, .true., Unassociated volume) + else + NN = NAdof(vvol) + + SALLOCATE( solution, (1:NN, 0:2), zero) + + ! Pack field and its derivatives + packorunpack = 'P' + WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,0), 0 ) ) ! packing; + WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,2), 2 ) ) ! packing; + + ! compute the field with renewed dpflux via single Newton method step + solution(1:NN, 0) = solution(1:NN, 0) - dpfluxout(Mvol) * solution(1:NN, 2) + + ! Unpack field in vector potential Fourier harmonics + packorunpack = 'U' + WCALL( dforce, packab, ( packorunpack, vvol, NN, solution(1:NN,0), 0 ) ) ! unpacking; + + DALLOCATE( solution ) + endif + endif + DALLOCATE(Fvec) DALLOCATE(dpfluxout) @@ -393,6 +439,28 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) lBBintegral(1:Nvol) = lBBintegral(1:Nvol) * half + voltotal = 0.0 + + do vvol = 1, Mvol + WCALL(dforce, IsMyVolume, (vvol)) + if( IsMyVolumeValue .EQ. 0 ) then + cycle + else if( IsMyVolumeValue .EQ. -1) then + FATAL(dforce, .true., Unassociated volume) + endif + + vvolume(vvol) = vvolume(vvol) + + press = adiabatic(vvol) * pscale / vvolume(vvol)**gamma + + betavol(vvol) = press * vvolume(vvol) / lBBintegral(vvol) + betavol(vvol) = betavol(vvol) * vvolume(vvol) + voltotal = voltotal+vvolume(vvol) + !write(*,*) "Calc beta: ", betavol(vvol), vvolume(vvol) + enddo + + BetaTotal = sum(betavol(1:Nvol))/voltotal ! Calculate total beta which is obtained from individual betas + Energy = sum( lBBintegral(1:Nvol) ) ! should also compute beta; !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! @@ -901,6 +969,17 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis) !call MPI_BARRIER( MPI_COMM_WORLD, ierr ) !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + + if(Lhessianallocated .and. Igeometry.eq.1) then + if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ; + write(*,*) "Writing .hessian file..." + open(munit, file=trim(ext)//".sp.hessian", status="unknown", form="unformatted") + write(munit) NGdof + write(munit) hessian(1:NGdof,1:NGdof) + close(munit) + + endif + endif RETURN(dforce) diff --git a/src/dfp100.f90 b/src/dfp100.f90 index acb75e37..2abd1b8a 100644 --- a/src/dfp100.f90 +++ b/src/dfp100.f90 @@ -63,7 +63,7 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives) DDtzcc, DDtzcs, DDtzsc, DDtzss, & DDzzcc, DDzzcs, DDzzsc, DDzzss, & dMA, dMB, dMD, dMG, MBpsi, solution, & - Nt, Nz, LILUprecond, Lsavedguvij, NOTMatrixFree, guvijsave, izbs + Nt, Nz, LILUprecond, Lsavedguvij, NOTMatrixFree, guvijsave, izbs, total_pflux LOCALS !------ @@ -205,6 +205,19 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives) Fvec(1:Mvol-1) = IPDt(1:Mvol-1) - Isurf(1:Mvol-1) endif + ! Set the total pflux to 0 + if(Igeometry.eq.1) then + + IPDtDpf(1,Mvol) = -pi2 * Bt00(1,1,2) + IPdtdPf(2:Mvol,Mvol) = 0.0 + IPDtDpf(Mvol,1:Mvol) = 1.0 + + ! Constraint the total pflux (pflux(Mvol)) + ! zero for symmetric initial profiles + ! non-zero for asymmetric profiles + Fvec(Mvol) = sum(dpflux(1:Mvol)) - total_pflux + endif + ! Compute poloidal linking current constraint as well in case of free boundary computation if ( Lfreebound.eq.1 ) then diff --git a/src/global.f90 b/src/global.f90 index b02862b1..bc5ae037 100644 --- a/src/global.f90 +++ b/src/global.f90 @@ -251,11 +251,14 @@ module allglobal REAL :: ForceErr !< total force-imbalance REAL :: Energy !< MHD energy REAL :: BnsErr !< (in freeboundary) error in self-consistency of field on plasma boundary (Picard iteration) + REAL :: BetaTotal = 0.0 !< Beta, averaged over entire domain REAL , allocatable :: IPDt(:), IPDtDpf(:,:) !< Toroidal pressure-driven current INTEGER :: Mvol + REAL :: total_pflux ! used when Lconstraint=3, Igeometry=1 + LOGICAL :: YESstellsym !< internal shorthand copies of Istellsym, which is an integer input; LOGICAL :: NOTstellsym !< internal shorthand copies of Istellsym, which is an integer input; @@ -1675,7 +1678,7 @@ subroutine wrtend write(iunit,'(" adiabatic = ",257es23.15)') adiabatic(1:Mvol) write(iunit,'(" mu = ",257es23.15)') mu(1:Mvol) write(iunit,'(" Ivolume = ",257es23.15)') Ivolume(1:Mvol) - write(iunit,'(" Isurf = ",257es23.15)') Isurf(1:Mvol-1), 0.0 + write(iunit,'(" Isurf = ",257es23.15)') IPDt(1:Mvol) ! Prints the actual surf current, not the targeted one write(iunit,'(" Lconstraint = ",i9 )') Lconstraint write(iunit,'(" pl = ",257i23 )') pl(0:Mvol) write(iunit,'(" ql = ",257i23 )') ql(0:Mvol) diff --git a/src/hesian.f90 b/src/hesian.f90 index a9a15b81..0c821c4f 100644 --- a/src/hesian.f90 +++ b/src/hesian.f90 @@ -22,7 +22,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) use inputlist, only : Wmacros, Whesian, Igeometry, Nvol, pflux, helicity, mu, Lfreebound, & LHevalues, LHevectors, LHmatrix, & Lperturbed, dpp, dqq, & - Lcheck, Lfindzero + Lcheck, Lfindzero, Lconstraint use cputiming, only : Thesian @@ -37,7 +37,9 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) Lhessian3Dallocated,denergydrr, denergydrz,denergydzr,denergydzz, & LocalConstraint -!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! + use sphdf5, only : write_stability + + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! LOCALS @@ -51,7 +53,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) REAL :: xx(0:NGdof,-2:2), ff(0:NGdof,-2:2), df(1:NGdof)!, deriv INTEGER :: vvol, idof, ii, mi, ni, irz, issym, isymdiff, lvol, ieval(1:1), igdof, ifd - REAL :: oldEnergy(-2:2), error, cpul + REAL :: oldEnergy(-2:2), error, cpul REAL :: oldBB(1:Mvol,-2:2), oBBdRZ(1:Mvol,0:1,1:LGdof), ohessian(1:NGdof,1:NGdof) @@ -90,6 +92,11 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) BEGIN(hesian) + ! Only makes sense to compute the Hessian if helicity is constrained + if(Lconstraint.ne.2) then + return + endif + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! lmu(1:Nvol) = mu(1:Nvol) ; lpflux(1:Nvol) = pflux(1:Nvol) ; lhelicity(1:Nvol) = helicity(1:Nvol) ! save original profile information; 20 Jun 14; @@ -229,7 +236,7 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) !if (LHmatrix) then Lhessian3Dallocated = .true. !else - ! Lhessianallocated = .true. + Lhessianallocated = .true. !endif !This step cleared. @@ -585,6 +592,10 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof ) endif ! end of if( ( LHevalues .or. LHevectors ) ) + ! output hessian, eigenvalues, and eigenvectors to the .h5 file + if( LHmatrix .and. Lconstraint.eq.2) then + WCALL( hesian, write_stability, (ohessian, evalr, evali, evecr, NGdof) ) + endif !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! if( myid.eq.0 .and. Lperturbed.eq.1 ) then diff --git a/src/newton.f90 b/src/newton.f90 index ee2fadde..25655733 100644 --- a/src/newton.f90 +++ b/src/newton.f90 @@ -97,7 +97,7 @@ subroutine newton( NGdof, position, ihybrd ) INTEGER :: ML, MU ! required for only Lc05ndf; - LOGICAL :: Lexit = .true. ! perhaps this could be made user input; + LOGICAL :: Lexit = .false. ! if false, hessian is printed even with initial force balance LOGICAL :: LComputeAxis INTEGER :: nprint = 1, nfev, njev diff --git a/src/preset.f90 b/src/preset.f90 index db447429..02e4c86e 100644 --- a/src/preset.f90 +++ b/src/preset.f90 @@ -508,8 +508,14 @@ subroutine preset SALLOCATE( dtflux, (1:Mvol), zero ) SALLOCATE( dpflux, (1:Mvol), zero ) + ! with Igeom=1, Lcons=3, need to provide total pflux + if(Lconstraint.eq.3 .and. Igeometry.eq.1) then + total_pflux = pflux(Mvol) * phiedge / pi2 + pflux(Mvol) = 0 + endif + select case( Igeometry ) - case( 1 ) ; dtflux(1) = tflux(1) ; dpflux(1) = pflux(1) ! Cartesian ; this is the "inverse" operation defined in xspech; 09 Mar 17; + case( 1 ) ; dtflux(1) = tflux(1) ; dpflux(1) = pflux(1) ! Cartesian ; this is the "inverse" operation defined in xspech; 09 Mar 17; case( 2:3 ) ; dtflux(1) = tflux(1) ; dpflux(1) = zero ! cylindrical or toroidal; end select @@ -807,7 +813,12 @@ subroutine preset if( Lfreebound.eq.1 ) then SALLOCATE( IPDtDpf, (1:Mvol , 1:Mvol ), zero) else - SALLOCATE( IPDtDpf, (1:Mvol-1, 1:Mvol-1), zero) + if(Igeometry.eq.1) then + ! add an additional constraint to make the total pflux = 0 + SALLOCATE( IPDtDpf, (1:Mvol, 1:Mvol), zero) + else + SALLOCATE( IPDtDpf, (1:Mvol-1, 1:Mvol-1), zero) + endif endif SALLOCATE( cheby, (0:Mrad,0:2), zero ) diff --git a/src/sphdf5.f90 b/src/sphdf5.f90 index 756414c9..45575bf2 100644 --- a/src/sphdf5.f90 +++ b/src/sphdf5.f90 @@ -45,6 +45,7 @@ module sphdf5 integer(hid_t) :: dt_nDcalls_id !< Memory datatype identifier (for "nDcalls" dataset in "/iterations") integer(hid_t) :: dt_Energy_id !< Memory datatype identifier (for "Energy" dataset in "/iterations") integer(hid_t) :: dt_ForceErr_id !< Memory datatype identifier (for "ForceErr" dataset in "/iterations") + integer(hid_t) :: dt_BetaTotal_id !< Memory datatype identifier (for "ForceErr" dataset in "/iterations") integer(hid_t) :: dt_iRbc_id !< Memory datatype identifier (for "iRbc" dataset in "/iterations") integer(hid_t) :: dt_iZbs_id !< Memory datatype identifier (for "iZbs" dataset in "/iterations") integer(hid_t) :: dt_iRbs_id !< Memory datatype identifier (for "iRbs" dataset in "/iterations") @@ -480,6 +481,7 @@ subroutine init_convergence_output H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, type_size_i, dt_nDcalls_id, hdfier) ) H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, type_size_d, dt_Energy_id, hdfier) ) H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, type_size_d, dt_ForceErr_id, hdfier) ) + H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, type_size_d, dt_BetaTotal_id, hdfier) ) H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, irbc_size, dt_iRbc_id, hdfier) ) H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, irbc_size, dt_iZbs_id, hdfier) ) H5CALL( sphdf5, h5tcreate_f, (H5T_COMPOUND_F, irbc_size, dt_iRbs_id, hdfier) ) @@ -488,6 +490,7 @@ subroutine init_convergence_output H5CALL( sphdf5, h5tinsert_f, (dt_nDcalls_id, "nDcalls", offset, H5T_NATIVE_INTEGER, hdfier) ) H5CALL( sphdf5, h5tinsert_f, (dt_Energy_id, "Energy", offset, H5T_NATIVE_DOUBLE, hdfier) ) H5CALL( sphdf5, h5tinsert_f, (dt_ForceErr_id, "ForceErr", offset, H5T_NATIVE_DOUBLE, hdfier) ) + H5CALL( sphdf5, h5tinsert_f, (dt_BetaTotal_id, "BetaTotal", offset, H5T_NATIVE_DOUBLE, hdfier) ) H5CALL( sphdf5, h5tinsert_f, (dt_iRbc_id, "iRbc", offset, iRZbscArray_id, hdfier) ) H5CALL( sphdf5, h5tinsert_f, (dt_iZbs_id, "iZbs", offset, iRZbscArray_id, hdfier) ) H5CALL( sphdf5, h5tinsert_f, (dt_iRbs_id, "iRbs", offset, iRZbscArray_id, hdfier) ) @@ -511,7 +514,7 @@ end subroutine init_convergence_output !> subroutine write_convergence_output( nDcalls, ForceErr ) - use allglobal, only : myid, mn, Mvol, Energy, iRbc, iZbs, iRbs, iZbc + use allglobal, only : myid, mn, Mvol, Energy, iRbc, iZbs, iRbs, iZbc, BetaTotal LOCALS INTEGER, intent(in) :: nDcalls @@ -545,6 +548,8 @@ subroutine write_convergence_output( nDcalls, ForceErr ) & mem_space_id=memspace, file_space_id=dataspace, xfer_prp=plist_id), __FILE__, __LINE__) H5CALL( sphdf5, h5dwrite_f, (iteration_dset_id, dt_ForceErr_id, ForceErr, INT((/1/), HSIZE_T), hdfier, & & mem_space_id=memspace, file_space_id=dataspace, xfer_prp=plist_id), __FILE__, __LINE__) + H5CALL( sphdf5, h5dwrite_f, (iteration_dset_id, dt_BetaTotal_id, BetaTotal, INT((/1/), HSIZE_T), hdfier, & + & mem_space_id=memspace, file_space_id=dataspace, xfer_prp=plist_id), __FILE__, __LINE__) H5CALL( sphdf5, h5dwrite_f, (iteration_dset_id, dt_iRbc_id, iRbc, INT((/mn,Mvol+1/), HSIZE_T), hdfier, & & mem_space_id=memspace, file_space_id=dataspace, xfer_prp=plist_id), __FILE__, __LINE__) H5CALL( sphdf5, h5dwrite_f, (iteration_dset_id, dt_iZbs_id, iZbs, INT((/mn,Mvol+1/), HSIZE_T), hdfier, & @@ -977,12 +982,50 @@ end subroutine write_vector_potential !> \brief Write the final state of the equilibrium to the output file. !> \ingroup grp_output !> + +subroutine write_stability(ohessian, evalr, evali, evecr, NGdof) + + use inputlist, only : LHevalues, LHevectors, LHmatrix + + LOCALS + integer, intent(in) :: NGdof + REAL, intent(in) :: ohessian(:,:), evalr(:), evali(:), evecr(:,:) + integer(hid_t) :: grpStability + + BEGIN( sphdf5 ) + + if (.not.LHevalues .and. .not.LHevectors .and. .not.LHmatrix) return + + if (myid.eq.0 .and. .not.skip_write) then + + HDEFGRP( file_id, stability, grpStability) + + if (LHmatrix) then + HWRITERA( grpStability, NGdof, NGdof, ohessian, ohessian(1:NGdof,1:NGdof) ) + endif + + if (LHevectors) then + HWRITERA( grpStability, NGdof, NGdof, evecr, evecr(1:NGdof,1:NGdof) ) + endif + + if (LHevalues) then + HWRITERV( grpStability, NGdof, evalr, evalr(1:NGdof) ) + HWRITERV( grpStability, NGdof, evali, evali(1:NGdof) ) + endif + + HCLOSEGRP( grpStability ) + + endif ! myid.eq.0 + +end subroutine write_stability + + subroutine hdfint use fileunits, only : ounit use inputlist use allglobal, only : ncpu, cpus, & - Mvol, ForceErr, BnsErr,& + Mvol, ForceErr, BnsErr, BetaTotal, & mn, im, in, iRbc, iZbs, iRbs, iZbc, & mns, ims, ins, & dRbc, dZbs, dRbs, dZbc, & @@ -1046,6 +1089,8 @@ subroutine hdfint HWRITERV( grpOutput, 1, BnsErr, (/ BnsErr /)) ! already in /input/global !latex \type{ForceErr} & real & \pb{force-balance error across interfaces} \\ HWRITERV( grpOutput, 1, ForceErr, (/ ForceErr /)) +!latex \type{BetaTotal} & real & \pb{Total plasma beta} \\ + HWRITERV( grpOutput, 1, BetaTotal, (/ BetaTotal /)) !latex \type{Ivolume} & real & \pb{Volume current at output (parallel, externally induced)} HWRITERV( grpOutput, Mvol, Ivolume, Ivolume(1:Mvol)) !latex \type{IPDt} & real & \pb{Surface current at output} @@ -1150,6 +1195,7 @@ subroutine finish_outfile H5CALL( sphdf5, h5tclose_f, (dt_nDcalls_id, hdfier) , __FILE__, __LINE__) H5CALL( sphdf5, h5tclose_f, (dt_Energy_id, hdfier) , __FILE__, __LINE__) H5CALL( sphdf5, h5tclose_f, (dt_ForceErr_id, hdfier) , __FILE__, __LINE__) + H5CALL( sphdf5, h5tclose_f, (dt_BetaTotal_id, hdfier) , __FILE__, __LINE__) H5CALL( sphdf5, h5tclose_f, (dt_iRbc_id, hdfier) , __FILE__, __LINE__) H5CALL( sphdf5, h5tclose_f, (dt_iZbs_id, hdfier) , __FILE__, __LINE__) H5CALL( sphdf5, h5tclose_f, (dt_iRbs_id, hdfier) , __FILE__, __LINE__) diff --git a/src/volume.f90 b/src/volume.f90 index 76581dd1..5ee2e635 100644 --- a/src/volume.f90 +++ b/src/volume.f90 @@ -41,7 +41,7 @@ subroutine volume( lvol, vflag ) use fileunits, only : ounit - use inputlist, only : Wvolume, Igeometry, Nvol, pscale + use inputlist, only : Wvolume, Igeometry, Nvol, pscale, rpol, rtor use cputiming @@ -120,6 +120,9 @@ subroutine volume( lvol, vflag ) if( dBdX%issym.eq.0 ) dvolume = one ! note that the sign factor for the lower interface is included below; 20 Jun 14; endif + ! respecting rpol, rtor sizes of slab + vol(innout) = vol(innout) * rpol * rtor + !-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-! case( 2 ) !>
  • \c Igeometry.eq.2 : cylindrical : \f$\sqrt g = R R_s = \frac{1}{2}\partial_s (R^2)\f$