Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

these modifications compute the virtual-casing self-consistent vacuum… #135

Open
wants to merge 3 commits into
base: master
Choose a base branch
from
Open

these modifications compute the virtual-casing self-consistent vacuum… #135

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
4de43fc
these modifications compute the virtual-casing self-consistent vacuum…
SRHudson Apr 9, 2021
56f68a9
laplace-virtual-casing matrix implemented
SRHudson Jun 29, 2022
b4995e0
Implemented derivatives of dVc; work in progress
abaillod Aug 23, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 7 additions & 5 deletions bfield.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ subroutine bfield( zeta, st, Bst ) ! the format of this subroutine is constraine

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

use constants, only : zero, one, half, two
use constants, only : zero, one, half, two, ten

use numerical, only : vsmall, small

Expand Down Expand Up @@ -145,12 +145,14 @@ subroutine bfield( zeta, st, Bst ) ! the format of this subroutine is constraine

if( abs(gBzeta).lt.vsmall ) then

cput = GETTIME
! cput = GETTIME

write(ounit,'("bfield : ",f10.2," : lvol=",i3," ; zeta="es23.15" ; (s,t)=("es23.15" ,"es23.15" ) ; B^z="es23.15" ;")') &
cput-cpus, lvol, zeta, st(1:2), dBu(3)
! write(ounit,'("bfield : ",f10.2," : lvol=",i3," ; zeta="es23.15" ; (s,t)=("es23.15" ,"es23.15" ) ; B^z="es23.15" ;")') &
! cput-cpus, lvol, zeta, st(1:2), dBu(3)

! FATAL( bfield, abs(dBu(3)).lt.vsmall, field is not toroidal )

FATAL( bfield, abs(dBu(3)).lt.vsmall, field is not toroidal )
gBzeta = ten

endif

Expand Down
4 changes: 4 additions & 0 deletions df00ab.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,8 @@ subroutine df00ab( pNN , xi , Fxi , DFxi , Ldfjac , iflag )
BEGIN(df00ab)

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

! write(ounit,'("df00ab : ", 10x ," : 1.1000 xi =",99(es23.15,","))') xi(0)

#ifdef DEBUG

Expand Down Expand Up @@ -80,6 +82,8 @@ subroutine df00ab( pNN , xi , Fxi , DFxi , Ldfjac , iflag )

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

! write(ounit,'("df00ab : ", 10x ," : 1.9000 xi =",99(es23.15,","))') xi(0)

RETURN(df00ab)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Expand Down
50 changes: 44 additions & 6 deletions dforce.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,7 +107,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
Lconstraint, Lcheck, dRZ, &
Lextrap, &
mupftol, &
Lfreebound
Lfreebound, &
mu

use cputiming, only : Tdforce

Expand Down Expand Up @@ -172,6 +173,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

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

! write(ounit,'("dforce : ", 10x ," : 2.2000 mu =",99(es23.15,","))') mu(1:Nvol)

! Unpack position to generate arrays iRbc, iZbs, IRbs, iZbc.

packorunpack = 'U' ! unpack geometrical degrees-of-freedom;
Expand All @@ -183,6 +186,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
WCALL( dforce, packxi,( NGdof, position(0:NGdof), Mvol, mn, iRbc(1:mn,0:Mvol), iZbs(1:mn,0:Mvol), &
iRbs(1:mn,0:Mvol), iZbc(1:mn,0:Mvol), packorunpack, LcomputeDerivatives, LComputeAxis ) )

! write(ounit,'("dforce : ", 10x ," : 2.2100 mu =",99(es23.15,","))') mu(1:Nvol)

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

if( LcomputeDerivatives ) then
Expand All @@ -192,7 +197,6 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
dBBdmp(1:LGdof,1:Mvol,0:1,1:2) = zero
endif


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

! SOLVE FIELD IN AGREEMENT WITH CONSTRAINTS AND GEOMETRY
Expand All @@ -218,15 +222,23 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

! Solve for field
dBdX%L = LComputeDerivatives

! write(ounit,'("dforce : ", 10x ," : 2.2110 mu =",99(es23.15,","))') mu(1:Nvol)

WCALL(dforce, dfp100, (Ndofgl, Xdof, Fvec, iflag) )

! write(ounit,'("dforce : ", 10x ," : 2.2120 mu =",99(es23.15,","))') mu(1:Nvol)

DALLOCATE( Fvec )

!do vvol=1,Mvol
! WCALL(dforce, brcast, (vvol) )
!enddo
! --------------------------------------------------------------------------------------------------
! Global constraint - call the master thread calls hybrd1 on dfp100, others call dfp100_loop.

! write(ounit,'("dforce : ", 10x ," : 2.2200 mu =",99(es23.15,","))') mu(1:Nvol)

else

IPDtdPf = zero
Expand All @@ -242,14 +254,20 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

SALLOCATE( Fvec, (1:Ndofgl), zero )

! write(ounit,'("dforce : ", 10x ," : 2.2210 mu =",99(es23.15,","))') mu(1:Nvol)

WCALL(dforce, dfp100, (Ndofgl, Xdof(1:Mvol-1), Fvec(1:Ndofgl), 1))

SALLOCATE(dpfluxout, (1:Ndofgl), zero )
if ( myid .eq. 0 ) then

! write(ounit,'("dforce : ", 10x ," : 2.2220 mu =",99(es23.15,","))') mu(1:Nvol)

dpfluxout = Fvec
call DGESV( Ndofgl, 1, IPdtdPf, Ndofgl, ipiv, dpfluxout, Ndofgl, idgesv )

! write(ounit,'("dforce : ", 10x ," : 2.2230 mu =",99(es23.15,","))') mu(1:Nvol)

! one step Newton's method
dpflux(2:Mvol) = dpflux(2:Mvol) - dpfluxout(1:Mvol-1)
if( Lfreebound.eq.1 ) then
Expand Down Expand Up @@ -318,8 +336,12 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
! end select
! #endif

! write(ounit,'("dforce : ", 10x ," : 2.2300 mu =",99(es23.15,","))') mu(1:Nvol)

endif !matches if( LocalConstraint )

! write(ounit,'("dforce : ", 10x ," : 2.2400 mu =",99(es23.15,","))') mu(1:Nvol)

! --------------------------------------------------------------------------------------------------
! MPI COMMUNICATIONS

Expand Down Expand Up @@ -354,12 +376,18 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
enddo

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

! write(ounit,'("dforce : ", 10x ," : 2.2400 mu =",99(es23.15,","))') mu(1:Nvol)

! Compute local force and derivatives
WCALL(dforce, dfp200, ( LcomputeDerivatives, vvol) )

! write(ounit,'("dforce : ", 10x ," : 2.2500 mu =",99(es23.15,","))') mu(1:Nvol)

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

#ifdef DEBUG

if( Lcheck.eq.2 ) then
write(ounit,'("dforce : ", 10x ," : myid=",i3," ; finished computing derivatives of rotational-transform wrt mu and dpflux ;")') myid
stop "dforce : : myid= ; finished computing derivatives of rotational-transform wrt mu and dpflux ;" ! this will allow other cpus to finish;
Expand All @@ -368,8 +396,11 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)
FATAL( dforce, Lcheck.eq.2, finished computing derivatives of rotational-transform wrt mu and dpflux )

if( Wdforce ) write(ounit,'("dforce : " 10x " : myid="i3" ; LComputeDerivatives="L2" ; ImagneticOK="999L2)') myid, LComputeDerivatives, ImagneticOK(1:Mvol)

#endif

! write(ounit,'("dforce : ", 10x ," : 2.2500 mu =",99(es23.15,","))') mu(1:Nvol)

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

! Broadcast information to all CPUs
Expand All @@ -394,6 +425,7 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

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

! write(ounit,'("dforce : ", 10x ," : 2.2600 mu =",99(es23.15,","))') mu(1:Nvol)

! CONSTRUCT FORCE
! ---------------
Expand Down Expand Up @@ -481,6 +513,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

enddo ! end of do vvol;

! write(ounit,'("dforce : ", 10x ," : 2.2700 mu =",99(es23.15,","))') mu(1:Nvol)

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

if( NGdof.ne.0 ) then ; ForceErr = sqrt( sum( force(1:NGdof)*force(1:NGdof) ) / NGdof ) ! this includes spectral constraints;
Expand All @@ -505,6 +539,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

#endif

! write(ounit,'("dforce : ", 10x ," : 2.2800 mu =",99(es23.15,","))') mu(1:Nvol)

4000 format("dforce : ",f10.2," : ",6x,3x,"; ",:,"|f|=",es12.5," ; ",:,"time=",f10.2,"s ;",:," log",a5,"=",28f6.2 ," ...")
4001 format("dforce : ", 10x ," : ",6x,3x,"; ",:," ", 12x ," ",:," ", 10x ," ;",:," log",a5,"=",28f6.2 ," ...")

Expand Down Expand Up @@ -678,6 +714,8 @@ subroutine dforce( NGdof, position, force, LComputeDerivatives, LComputeAxis)

!call MPI_BARRIER( MPI_COMM_WORLD, ierr )

! write(ounit,'("dforce : ", 10x ," : 2.2900 mu =",99(es23.15,","))') mu(1:Nvol)

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

RETURN(dforce)
Expand Down Expand Up @@ -817,7 +855,7 @@ subroutine fndiff_dforce( NGdof )
write(ounit,'(A)') NEW_LINE('A')

do ii=1, NGdof
write(ounit,1345) myid, im(ii), in(ii), hessian(ii,:)
write(ounit,1345) myid, hessian(ii,:)
write(10 ,1347) hessian(ii,:)
enddo
close(10)
Expand All @@ -827,14 +865,14 @@ subroutine fndiff_dforce( NGdof )
! Print finite differences
open(10, file=trim(ext)//'.Lcheck6_output.FiniteDiff.txt', status='unknown')
do ii=1, NGdof
write(ounit,1347) myid, im(ii), in(ii), finitediff_estimate(ii,:)
write(ounit,1346) myid, finitediff_estimate(ii,:)
write(10 ,1347) finitediff_estimate(ii,:)
enddo
write(ounit,'(A)') NEW_LINE('A')
close(10)

1345 format("dforce: myid=",i3," ; (",i4,",",i4," ; Hessian = ",512f16.10 " ;")
1346 format("dforce: myid=",i3," ; (",i4,",",i4," ; Finite differences = ",512f16.10 " ;")
1345 format("dforce: myid=",i3,"; Hessian = ",512f16.10 " ;")
1346 format("dforce: myid=",i3,"; Finite differences = ",512f16.10 " ;")
1347 format(512F22.16, " ")

endif
Expand Down
17 changes: 11 additions & 6 deletions dfp100.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,16 +30,14 @@
!latex \subsection{Evaluation of constraint derivatives}
!latex Not implemented for now.



subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)

use constants, only : zero, half, one, two, pi2, pi, mu0

use fileunits, only : ounit

use inputlist, only : Wmacros, Wdfp100, Igeometry, Nvol, Lrad, Isurf, &
Lconstraint, Lfreebound, curpol
Lconstraint, Lfreebound, curpol, mu

use cputiming, only : Tdfp100

Expand Down Expand Up @@ -77,18 +75,18 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
LOGICAL :: LComputeDerivatives
INTEGER :: deriv, Lcurvature



BEGIN(dfp100)

! write(ounit,'("dfp100 : ", 10x ," : 2.2110 mu =",99(es23.15,","))') mu(1:Nvol)

dpflux(2:Mvol) = x - xoffset

! Now each CPU perform the calculation in its volume(s)
do vvol = 1, Mvol

LREGION(vvol) ! assigns Lcoordinatesingularity, Lplasmaregion, etc. ;
ImagneticOK(vvol) = .false.
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

! Determines if this volume vvol should be computed by this thread.
call IsMyVolume(vvol)
Expand Down Expand Up @@ -126,7 +124,11 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
WCALL( dfp100, matrixBG, ( vvol, mn, ll ) )
endif

! write(ounit,'("dfp100 : ", 10x ," : 2.2120 mu =",99(es23.15,","))') mu(1:Nvol)

WCALL( dfp100, ma02aa, ( vvol, NN ) )

! write(ounit,'("dfp100 : ", 10x ," : 2.2130 mu =",99(es23.15,","))') mu(1:Nvol)

Lsavedguvij = .false.

Expand Down Expand Up @@ -164,6 +166,7 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
endif
enddo

! write(ounit,'("dfp100 : ", 10x ," : 2.2140 mu =",99(es23.15,","))') mu(1:Nvol)

! Evaluation of global constraint and communications
if( .not.LocalConstraint ) then
Expand Down Expand Up @@ -223,6 +226,8 @@ subroutine dfp100(Ndofgl, x, Fvec, LComputeDerivatives)
end select
endif

! write(ounit,'("dfp100 : ", 10x ," : 2.2190 mu =",99(es23.15,","))') mu(1:Nvol)

RETURN(dfp100)

end subroutine dfp100
Loading