Skip to content

Commit

Permalink
added updates
Browse files Browse the repository at this point in the history
  • Loading branch information
@ErolBa
ErolBa committed Dec 12, 2024
1 parent 4f1d0da commit e1183cc
Show file tree
Hide file tree
Showing 8 changed files with 200 additions and 35 deletions.
89 changes: 84 additions & 5 deletions src/dforce.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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 )
Expand All @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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;

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Expand Down Expand Up @@ -902,6 +970,17 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

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)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Expand Down
15 changes: 14 additions & 1 deletion src/dfp100.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
!------
Expand Down Expand Up @@ -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

Expand Down
5 changes: 4 additions & 1 deletion src/global.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 !< total number of volumes (including the vacuum region in the case of free-boundary calculations)

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;

Expand Down Expand Up @@ -1678,7 +1681,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)
Expand Down
54 changes: 32 additions & 22 deletions src/hesian.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -37,6 +37,8 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )
Lhessian3Dallocated,denergydrr, denergydrz,denergydzr,denergydzz, &
LocalConstraint

use sphdf5, only : write_stability

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

LOCALS
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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.

Expand Down Expand Up @@ -416,16 +423,16 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

if( LHmatrix ) then

if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ;
open(munit, file=trim(get_hidden(ext))//".GF.ma", status="unknown", form="unformatted")
write(munit) NGdof
write(munit) ohessian(1:NGdof,1:NGdof)
close(munit)
endif

endif
! if( LHmatrix ) then
!
! if( myid.eq.0 ) then ; cput = GETTIME ; write(ounit,'("hesian : ",f10.2," : LHmatrix="L2" ;")')cput-cpus, LHmatrix ;
! open(munit, file=trim(get_hidden(ext))//".GF.ma", status="unknown", form="unformatted")
! write(munit) NGdof
! write(munit) ohessian(1:NGdof,1:NGdof)
! close(munit)
! endif
!
! endif


! if( myid.eq.0 .and. ( LHevalues .or. LHevectors ) ) then ! the call to dforce below requires all cpus; 04 Dec 14;
Expand Down Expand Up @@ -572,19 +579,22 @@ subroutine hesian( NGdof, position, Mvol, mn, LGdof )

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!


if( myid.eq.0 ) then ! write to file; 04 Dec 14;
open(hunit, file=trim(get_hidden(ext))//".GF.ev", status="unknown", form="unformatted")
write(hunit) NGdof, Ldvr, Ldvi
write(hunit) evalr
write(hunit) evali
write(hunit) evecr
write(hunit) eveci
close(hunit)
endif ! end of if( myid.eq.0 ) ; 04 Dec 14;
! if( myid.eq.0 ) then ! write to file; 04 Dec 14;
! open(hunit, file=trim(get_hidden(ext))//".GF.ev", status="unknown", form="unformatted")
! write(hunit) NGdof, Ldvr, Ldvi
! write(hunit) evalr
! write(hunit) evali
! write(hunit) evecr
! write(hunit) eveci
! close(hunit)
! endif ! end of if( myid.eq.0 ) ; 04 Dec 14;

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
Expand Down
2 changes: 1 addition & 1 deletion src/newton.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
15 changes: 13 additions & 2 deletions src/preset.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -507,9 +507,15 @@ 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

Expand Down Expand Up @@ -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 )
Expand Down
Loading

0 comments on commit e1183cc

Please sign in to comment.