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.

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

Grid virtual casing #217

Open
wants to merge 31 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
31 commits
Select commit Hold shift + click to select a range
689370a
grid based virtual casing implemented
missing-user Nov 29, 2024
8339136
documentation
missing-user Nov 29, 2024
5d26ac4
accuracy estimates of virtual casing method
missing-user Nov 29, 2024
5502bb3
switch between grid based VC or adaptive VC
missing-user Nov 29, 2024
9ece12b
variable renaming, added to outputlist
missing-user Nov 30, 2024
e89bbcf
add to inputlist copy
missing-user Nov 30, 2024
87b10c4
unused variables
missing-user Nov 30, 2024
ce379ca
variable naming and documentation
missing-user Dec 7, 2024
f3a3fa4
Revert broken WIP code
missing-user Dec 7, 2024
20d9a04
rename to vcNt, vcNz
missing-user Dec 7, 2024
745e4c4
Revert derivative comment deletion
missing-user Dec 7, 2024
57a2d3f
typo prevented compilation
missing-user Dec 7, 2024
b2a3f0d
MPI parallel casinggrid
missing-user Dec 12, 2024
1e8ec20
typo in loop counter
missing-user Dec 12, 2024
544cdd0
logging only on rank0
missing-user Dec 12, 2024
497f187
Merge branch 'PrincetonUniversity:master' into grid-virtual-casing2
missing-user Dec 12, 2024
389b4b1
Exploit symmetries in VirtualCasing
missing-user Dec 12, 2024
dfff9a6
remove NaN check
missing-user Dec 12, 2024
f7adb3c
Fix casinggrid bug
missing-user Dec 13, 2024
e2c430b
documented, benchmarked, increased accuracy casinggrid
missing-user Dec 13, 2024
3e13472
Merge branch 'grid-virtual-casing2' of https://github.com/missing-use…
missing-user Dec 13, 2024
8c80363
correct stride regardless of Nt>Nz or Nt<Nz
missing-user Dec 13, 2024
7491933
max relative error instead of mean
missing-user Dec 14, 2024
e3cb69e
add to tests
missing-user Dec 16, 2024
461ebec
increased testcase resolution, improved logging
missing-user Dec 16, 2024
d8ecd62
missing zero initialization broke MPI, casinggrid can now skip parts …
missing-user Dec 16, 2024
78ff6dc
wrong build.yml xspec path
missing-user Dec 16, 2024
df4b737
reduced required resolution for gridVC
missing-user Dec 16, 2024
90cf04d
Merge branch 'grid-virtual-casing2' of https://github.com/missing-use…
missing-user Dec 16, 2024
0b6ab77
captalize vcNt vcNz
missing-user Dec 17, 2024
2cf26c2
update documentation
missing-user Dec 19, 2024
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
Prev Previous commit
Next Next commit
variable renaming, added to outputlist
  • Loading branch information
missing-user committed Dec 4, 2024
commit 9ece12bc7bb7dfbe57db19c9346e0c807a27e47c
30 changes: 11 additions & 19 deletions src/bnorml.f90
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )

use fileunits, only : ounit, lunit

use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad, lgridvcasing
use inputlist, only : Wmacros, Wbnorml, Igeometry, Lcheck, vcasingtol, vcasingper, Lrad, Lvcgrid

use cputiming, only : Tbnorml

Expand All @@ -99,7 +99,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
Nt, Nz, cfmn, sfmn, &
ijreal, ijimag, jireal, jiimag, &
globaljk, virtualcasingfactor, gteta, gzeta, Dxyz, Nxyz, &
Jxyz, Pbxyz, vcNtz
Jxyz, Pbxyz

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

Expand All @@ -114,6 +114,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
REAL :: accuracyestimate, resulth, resulth2, resulth4, deltah4h2, deltah2h
INTEGER :: vcstep


BEGIN(bnorml)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
Expand Down Expand Up @@ -144,14 +145,14 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
! When comparing results, both methods should always run
if (.true.) then
#else
if ( lgridvcasing.eq.1 ) then
if ( Lvcgrid.eq.1 ) then
#endif
! Precompute Jxyz(1:Ntz,1:3) and the corresponding positions on the high resolution plasma boundary

!$OMP PARALLEL DO SHARED(Pbxyz, Jxyz) PRIVATE(jk, teta, zeta)
!$OMP PARALLEL DO SHARED(Pbxyz, Jxyz) PRIVATE(jk, teta, zeta) COLLAPSE(2)
do kk = 0, vcNz-1 ;
zeta = kk * pi2nfp / vcNz
do jj = 0, vcNt-1 ;
zeta = kk * pi2nfp / vcNz
teta = jj * pi2 / vcNt ;
jk = 1 + jj + kk*vcNt

Expand All @@ -164,10 +165,10 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
do vcstep = 3, 0, -1
!$OMP PARALLEL DO SHARED(Dxyz, Pbxyz, Jxyz, ijimag) PRIVATE(jk, gBn)
do jk = 1, Ntz
call casing2( Dxyz(:,jk), Nxyz(:,jk), Pbxyz(1:vcNtz, 1:3), Jxyz(1:vcNtz, 1:3), 2**vcstep, gBn)
call casing2( Dxyz(:,jk), Nxyz(:,jk), Pbxyz, Jxyz, 2**vcstep, gBn)

ijreal(jk) = ijimag(jk)
ijimag(jk) = gBn
ijreal(jk) = ijimag(jk) ! previous solution (lower resolution)
ijimag(jk) = gBn ! current solution (higher resolution)

enddo
deltah4h2 = deltah2h
Expand All @@ -188,7 +189,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
! When comparing results, both methods should always run
if (.true.) then
#else
else ! if not lgridvcasing
else ! if not Lvcgrid
#endif

do kk = 0, Nz-1 ;
Expand Down Expand Up @@ -219,7 +220,7 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )

enddo ! end of do jj;
enddo ! end of do kk;
endif ! end of if (lgridvcasing )
endif ! end of if (Lvcgrid )
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

1001 format("bnorml : ", 10x ," : "a1" : (t,z) = ("f8.4","f8.4" ) ; gBn=",f23.15," ; ":" error =",f23.15" ;")
Expand All @@ -241,11 +242,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
case( 0 ) ! Lparallel = 0 ; 09 Mar 17;

RlBCAST(ijreal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp) ! plasma; 03 Apr 13;
!RlBCAST(ijimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)

!RlBCAST(jireal(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)
!RlBCAST(jiimag(1+kk*Nt:Nt+kk*Nt),Nt,kkmodnp)

case( 1 ) ! Lparallel = 1 ; 09 Mar 17;

do jj = 0, Nt-1
Expand All @@ -255,10 +251,6 @@ subroutine bnorml( mn, Ntz, efmn, ofmn )
jkmodnp = modulo(jk-1,ncpu)

RlBCAST(ijreal(jk),1,jkmodnp) ! plasma; 03 Apr 13;
!RlBCAST(ijimag(jk),1,jkmodnp)

!RlBCAST(jireal(jk),1,jkmodnp)
!RlBCAST(jiimag(jk),1,jkmodnp)

enddo

Expand Down
6 changes: 3 additions & 3 deletions src/casing.f90
Original file line number Diff line number Diff line change
Expand Up @@ -230,13 +230,13 @@ subroutine casing2( xyz, nxyz, Pbxyz, Jxyz, vcstep,gBn)

use inputlist, only : vcasingeps, vcNz, vcNt

use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC, vcNtz !, Pbxyz, Jxyz
use allglobal, only : myid, ncpu, cpus, MPI_COMM_SPEC !, Pbxyz, Jxyz

LOCALS

REAL, intent(in) :: xyz(3) ! arbitrary location; Cartesian;
REAL, intent(in) :: nxyz(3) ! surface normal on the computational boundary; Cartesian;
REAL, intent(in) :: Pbxyz(1:vcNtz, 1:3), Jxyz(1:vcNtz, 1:3)
REAL, intent(in) :: Pbxyz(1:vcNz*vcNt, 1:3), Jxyz(1:vcNz*vcNt, 1:3)
INTEGER, intent(in) :: vcstep
REAL, intent(out) :: gBn ! B.n on the computational boundary;

Expand All @@ -245,7 +245,7 @@ subroutine casing2( xyz, nxyz, Pbxyz, Jxyz, vcstep,gBn)

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
! plasmaNtz = SIZE(Pbxyz, 1)
plasmaNtz = vcNtz
plasmaNtz = vcNz*vcNt
! loop over the high resolution plasma boundary (inner boundary for virtual casing)
do jk = 1, plasmaNtz, vcstep ;
! position on computational boundary - position on plasma boundary
missing-user marked this conversation as resolved.
Show resolved Hide resolved
Expand Down
2 changes: 0 additions & 2 deletions src/global.f90
Original file line number Diff line number Diff line change
Expand Up @@ -844,8 +844,6 @@ module allglobal
INTEGER :: IBerror !< for computing error in magnetic field

INTEGER :: nfreeboundaryiterations !< number of free-boundary iterations already performed

INTEGER :: vcNtz !< number of grid points in \f$\theta\f$ and \f$\zeta\f$ for the plasma surface current computation
!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!

INTEGER, parameter :: Node = 2 !< best to make this global for consistency between calling and called routines
Expand Down
6 changes: 3 additions & 3 deletions src/inputlist.f90
Original file line number Diff line number Diff line change
Expand Up @@ -315,7 +315,7 @@ module inputlist
!< </ul>
INTEGER :: Ntoraxis = 3 !< the number of \f$n\f$ harmonics used in the Jacobian \f$m=1\f$ harmonic elimination method;
!< only relevant if \c Lrzaxis.ge.1 .
INTEGER :: Lgridvcasing = 0 !< Which method to use for the virtual casing integral. 0 = adaptive integration routine with guaranteed accuracy, 1 = fixed resolution grid
INTEGER :: Lvcgrid = 0 !< Which method to use for the virtual casing integral. 0 = adaptive integration routine with guaranteed accuracy, 1 = fixed resolution grid
!> @}

!> \addtogroup grp_global_local locallist
Expand Down Expand Up @@ -679,7 +679,7 @@ module inputlist
Mregular ,&
Lrzaxis ,&
Ntoraxis ,&
Lgridvcasing
Lvcgrid

namelist/locallist/&
LBeltrami ,&
Expand Down Expand Up @@ -895,7 +895,7 @@ subroutine initialize_inputs
Mregular = -1
Lrzaxis = 1
Ntoraxis = 3
Lgridvcasing = 0
Lvcgrid = 0

! locallist

Expand Down
7 changes: 3 additions & 4 deletions src/preset.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1723,17 +1723,16 @@ subroutine preset

!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!-!
if (Lfreebound > 0) then ! Only do for free-boundary; 7 Nov 18;
if (Lgridvcasing .eq. 1) then
if (Lvcgrid .eq. 1) then
if ( vcNt.lt.Nt ) then
FATAL( bnorml, .true., The plasma boundary resolution for virtual casing vcNt must be greater than or equal to Nt )
endif
if ( vcNz.lt.Nz ) then
FATAL( bnorml, .true., The plasma boundary resolution for virtual casing vcNz must be greater than or equal to Nz )
endif

vcNtz = vcNt*vcNz ! the plasma boundary has to use a higher resolution than the computational boundary
SALLOCATE( Jxyz, (1:vcNtz,1:3), zero ) ! Cartesian components of virtual casing surface current; needs to be recalculated at each iteration;
SALLOCATE( Pbxyz, (1:vcNtz,1:3), zero ) ! Cartesian points on the plasma boundary; needs to be recalculated at each iteration;
SALLOCATE( Jxyz, (1:vcNt*vcNz,1:3), zero ) ! Cartesian components of virtual casing surface current; needs to be recalculated at each iteration;
SALLOCATE( Pbxyz, (1:vcNt*vcNz,1:3), zero ) ! Cartesian points on the plasma boundary; needs to be recalculated at each iteration;
endif

SALLOCATE( Dxyz, (1:3,1:Ntz), zero ) ! Cartesian components of computational boundary; position; 14 Apr 17;
Expand Down
3 changes: 3 additions & 0 deletions src/sphdf5.f90
Original file line number Diff line number Diff line change
Expand Up @@ -294,6 +294,7 @@ subroutine mirror_input_to_outfile
HWRITEIV( grpInputNumerics, 1, Mregular , (/ Mregular /))
HWRITEIV( grpInputNumerics, 1, Lrzaxis , (/ Lrzaxis /))
HWRITEIV( grpInputNumerics, 1, Ntoraxis , (/ Ntoraxis /))
HWRITEIV( grpInputNumerics, 1, Lvcgrid , (/ Lvcgrid /))

HCLOSEGRP( grpInputNumerics, __FILE__, __LINE__)

Expand Down Expand Up @@ -344,6 +345,8 @@ subroutine mirror_input_to_outfile
HWRITEIV( grpInputGlobal, 1, vcasingits , (/ vcasingits /))
HWRITEIV( grpInputGlobal, 1, vcasingper , (/ vcasingper /))
HWRITEIV( grpInputGlobal, 1, mcasingcal , (/ mcasingcal /)) ! redundant;
HWRITEIV( grpInputGlobal, 1, vcnt , (/ vcnt /))
missing-user marked this conversation as resolved.
Show resolved Hide resolved
HWRITEIV( grpInputGlobal, 1, vcnz , (/ vcnz /))

HCLOSEGRP( grpInputGlobal )

Expand Down