BiotSavart.F90

!!matvec with matmul !#define __MATVEC_N(y,Mat,Vec) y=MATMUL(Mat,Vec) !#define __MATVEC_T(y,Mat,Vec) y=MATMUL(Vec,Mat) !#define __PMATVEC_N(fy,y,Mat,Vec) y=fyy+MATMUL(Mat,Vec) !#define __PMATVEC_T(fy,y,Mat,Vec) y=fyy+MATMUL(Vec,Mat) !#define __AMATVEC_N(y,fMat,Mat,Vec) y=fMatMATMUL(Mat,Vec) !#define __AMATVEC_T(y,fMat,Mat,Vec) y=fMatMATMUL(Vec,Mat) !#define __PAMATVEC_N(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Mat,Vec) !#define __PAMATVEC_T(fy,y,fMat,Mat,Vec) y=fyy+fMatMATMUL(Vec,Mat)

!!#define __GENERICMATVEC(NT,fy,y,fMat,Mat,Vec) CALL DGEMV(NT,SIZE(Mat,1),SIZE(Mat,2),fMat,Mat,SIZE(Mat,1),Vec,1,fy,y,1)

!!matmat with matmul !#define __MATMAT_NN(Y,A,B) Y=MATMUL(A,B) !#define __MATMAT_TN(Y,A,B) Y=MATMUL(TRANSPOSE(A),B) !#define __MATMAT_NT(Y,A,B) Y=MATMUL(A,TRANSPOSE(B)) !#define __MATMAT_TT(Y,A,B) Y=TRANSPOSE(MATMUL(B,A))

!#define __PMATMAT_NN(fy,Y,A,B) Y=fyY+MATMUL(A,B) !#define __PMATMAT_TN(fy,Y,A,B) Y=fyY+MATMUL(TRANSPOSE(A),B) !#define __PMATMAT_NT(fy,Y,A,B) Y=fyY+MATMUL(A,TRANSPOSE(B)) !#define __PMATMAT_TT(fy,Y,A,B) Y=fyY+TRANSPOSE(MATMUL(B,A))

!#define __AMATMAT_NN(Y,fa,A,B) Y=faMATMUL(A,B) !#define __AMATMAT_TN(Y,fa,A,B) Y=faMATMUL(TRANSPOSE(A),B) !#define __AMATMAT_NT(Y,fa,A,B) Y=faMATMUL(A,TRANSPOSE(B)) !#define __AMATMAT_TT(Y,fa,A,B) Y=faTRANSPOSE(MATMUL(B,A))

!#define __PAMATMAT_NN(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,B) !#define __PAMATMAT_TN(fy,Y,fa,A,B) Y=fyY+faMATMUL(TRANSPOSE(A),B) !#define __PAMATMAT_NT(fy,Y,fa,A,B) Y=fyY+faMATMUL(A,TRANSPOSE(B)) !#define __PAMATMAT_TT(fy,Y,fa,A,B) Y=fyY+faTRANSPOSE(MATMUL(B,A))

! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

!#define __GENERICMATMAT_NN(fy,Y,fa,A,B) CALL DGEMM('N','N',SIZE(A,1),SIZE(B,2),SIZE(B,1),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TN(fy,Y,fa,A,B) CALL DGEMM('T','N',SIZE(A,2),SIZE(B,2),SIZE(B,1),fa,A,SIZE(B,1),B,SIZE(B,1),fy,Y,SIZE(A,2)) !#define __GENERICMATMAT_NT(fy,Y,fa,A,B) CALL DGEMM('N','T',SIZE(A,1),SIZE(B,1),SIZE(B,2),fa,A,SIZE(A,1),B,SIZE(B,1),fy,Y,SIZE(A,1)) !#define __GENERICMATMAT_TT(fy,Y,fa,A,B) CALL DGEMM('T','T',SIZE(A,2),SIZE(B,1),SIZE(B,2),fa,A,SIZE(B,2),B,SIZE(B,1),fy,Y,SIZE(A,2))

! SIMPLE INTERFACE FOR DGEMM, specifying nrows/ncols of mat A and nrows/ncols of mat B (for any transpose!) ! GEMM does in general Y = fa A^?B^? + fy Y ! with structure: (m x n) = (m x k) (k x n) ! Y=A B : DGEMM('N','N',m,n,k,fa,Amat ,m, Bmat,k, fy,Y,m) ! Y=A^TB : DGEMM('T','N',m,n,k,fa,Amat ,k, Bmat,k, fy,Y,m) ! Y=A B^T : DGEMM('N','T',m,n,k,fa,Amat ,m, Bmat,n, fy,Y,m) ! Y=A^T*B^T : DGEMM('T','T',m,n,k,fa,Amat ,k, Bmat,n, fy,Y,m)

This file depends on

Help

Source Code

!===================================================================================================================================
! Copyright (c) 2025 GVEC Contributors, Max Planck Institute for Plasma Physics
! License: MIT
!===================================================================================================================================
#include "defines.FPP"

MODULE MODgvec_BiotSavart

IMPLICIT NONE
PUBLIC

CONTAINS

!================================================================================================================================!
!> Evaluate the magnetic field of a current loop discretized via line-segments.
!! The segments are defined by a list of `coil_points`, each segment being defined between points (i,i+1).
!! For a closed loop, the last point in the list is a repetition of the first point.
!! Each line segment is evaluated via the analytic compact Biot-Savart expression from
!! Hanson and Hirshman (2002) (https://doi.org/10.1063/1.1507589).
!! This implementation parallelizes over the number of evaluation positions `n_positions`.
!================================================================================================================================!
SUBROUTINE BiotSavart(n_positions, xyz, n_points, coil_points, prefactor, B)
    ! MODULES
    USE MODgvec_Globals, ONLY: TWOPI,CROSS,wp
    ! INPUT VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER,  INTENT(IN) :: n_positions               !! number of positions where the magnetic field is evaluated
    REAL(wp), INTENT(IN) :: xyz(3,n_positions)        !! x,y,z positions where the magnetic field is evaluated
    INTEGER,  INTENT(IN) :: n_points                  !! number of points representing the segmented coil
    REAL(wp), INTENT(IN) :: coil_points(3,n_points)   !! x,y,z positions of the segment chain (=polygon), n_segments=n_points-1
    REAL(wp), INTENT(IN) :: prefactor                 !! factor applied to the final magnetic field
    ! OUTPUT VARIABLES ----------------------------------------------------------------------------------------------------------!
    REAL(wp), INTENT(OUT):: B(3,n_positions)          !! magnetic field evaluated at xyz positions.
    ! LOCAL VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER              :: iPosition,iSegment,n_segments
    REAL(wp)             :: R_i(3)                    !! vector from current position to starting point of segment
    REAL(wp)             :: R_f(3)                    !! vector from current position to end point of segment (=starting point of next segment)
    REAL(wp)             :: mod_R_i,mod_R_f
    ! CODE ----------------------------------------------------------------------------------------------------------------------!
    n_segments = n_points-1
    B = 0.0_wp
    !$OMP PARALLEL DO     &
    !$OMP SCHEDULE(STATIC) DEFAULT(PRIVATE) SHARED(xyz, coil_points, n_positions, n_segments) &
    !$OMP REDUCTION(+:B)
    DO iPosition=1,n_positions
        R_f = xyz(:,iPosition)-coil_points(:,1)
        mod_R_f= SQRT(SUM(R_f*R_f))
        DO iSegment=1,n_segments
            R_i = R_f
            mod_R_i = mod_R_f
            R_f = xyz(:,iPosition)-coil_points(:,iSegment+1)
            mod_R_f = SQRT(SUM(R_f*R_f))
            !! RidotRf = SUM(R_i*R_f) ;  RixRf(:) = CROSS(R_i, R_f)
            !! segment_factor = (mod_R_i+mod_R_f)/(mod_R_i*mod_R_f*(mod_R_i*mod_R_f+RidotRf))
            !! B(:,iPosition) = B(:,iPosition)+segment_factor*RixRf
            B(:,iPosition) = B(:,iPosition)+((mod_R_i+mod_R_f)/(mod_R_i*mod_R_f*(mod_R_i*mod_R_f+SUM(R_i*R_f))))*CROSS(R_i, R_f)
        END DO !iSegment
    END DO !iPosition
    !$OMP END PARALLEL DO
    B = prefactor*B
END SUBROUTINE BiotSavart

!================================================================================================================================!
SUBROUTINE BiotSavart_VectorPotential(n_positions, xyz, n_points, n_segments, coil_points, ehat, L, prefactor, A)
    ! MODULES
    USE MODgvec_Globals, ONLY: TWOPI,CROSS,wp
    ! INPUT VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER,  INTENT(IN) :: n_positions               !! number of positions where the magnetic field is evaluated
    REAL(wp), INTENT(IN) :: xyz(3,n_positions)        !! x,y,z positions where the magnetic field is evaluated
    INTEGER,  INTENT(IN) :: n_points                  !! number of points representing the segmented coil
    INTEGER,  INTENT(IN) :: n_segments                !! number of coil segments
    REAL(wp), INTENT(IN) :: coil_points(3,n_points)   !! x,y,z positions of the segment chain (=polygon), n_segments=n_points-1
    REAL(wp), INTENT(IN) :: prefactor                 !! factor applied to the final magnetic field
    REAL(wp), INTENT(IN) :: ehat(3,n_segments)         !! unit vector along segment
    REAL(wp), INTENT(IN) :: L(n_segments)             !! segment length
    ! OUTPUT VARIABLES ----------------------------------------------------------------------------------------------------------!
    REAL(wp), INTENT(OUT):: A(3,n_positions)          !! magnetic field evaluated at xyz positions.
    ! LOCAL VARIABLES -----------------------------------------------------------------------------------------------------------!
    INTEGER              :: iPosition,iSegment
    REAL(wp)             :: R_i(3)                    !! vector from current position to starting point of segment
    REAL(wp)             :: R_f(3)                    !! vector from current position to end point of segment (=starting point of next segment)
    REAL(wp)             :: norm_ehat
    REAL(wp)             :: mod_R_i,mod_R_f
    REAL(wp)             :: f_epsilon
    ! CODE ----------------------------------------------------------------------------------------------------------------------!
    A = 0.0_wp
    !$OMP PARALLEL DO     &
    !$OMP SCHEDULE(STATIC) DEFAULT(PRIVATE) SHARED(xyz, coil_points, n_positions, n_segments, L, ehat) &
    !$OMP REDUCTION(+:A)
    DO iPosition=1,n_positions
        R_f = xyz(:,iPosition)-coil_points(:,1)
        mod_R_f= SQRT(SUM(R_f*R_f))
        DO iSegment=1,n_segments
            R_i = R_f
            mod_R_i = mod_R_f
            R_f = xyz(:,iPosition)-coil_points(:,iSegment+1)
            mod_R_f = SQRT(SUM(R_f*R_f))
            f_epsilon = LOG( (mod_R_i + mod_R_f + L(iSegment)) / (mod_R_i + mod_R_f - L(iSegment))  )
            A(:,iPosition) = A(:,iPosition) + ehat(:,iSegment) * f_epsilon
        END DO !iSegment
    END DO !iPosition
    !$OMP END PARALLEL DO
    A = prefactor*A
END SUBROUTINE BiotSavart_VectorPotential

END MODULE MODgvec_BiotSavart

BiotSavart.F90 Source File

Contents

Modules

Source Code

This file depends on

Source Code