diff --git a/src/fields_mod.f90 b/src/fields_mod.f90 index 817eef6..a2b0eec 100644 --- a/src/fields_mod.f90 +++ b/src/fields_mod.f90 @@ -1,1467 +1,1485 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: beam ! !> @author !> Patryk Kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for initializing the magnetic field, solving the Poisson equation and computing the moments of the particles distribution function !------------------------------------------------------------------------------ MODULE fields USE constants USE basic, ONLY: nr, nz, zgrid, rgrid, Br, Bz, Er, Ez, femorder, ngauss, nlppform, pot, Athet, & & splrz, splrz_ext, nlperiod, phinorm, nlPhis, nrank, mpirank, mpisize, step, it2d, timera, potxt, erxt, ezxt USE beam, ONLY: partslist USE bsplines USE mumps_bsplines use mpi Use omp_lib Use mpihelper, ONLY: db_type USE particletypes IMPLICIT NONE REAL(kind=db), allocatable, SAVE :: matcoef(:, :), phi_spline(:), vec1(:), vec2(:) REAL(kind=db), allocatable, SAVE :: loc_moments(:, :), loc_rhs(:), gradgtilde(:), fverif(:), ppformwork(:,:,:) INTEGER, SAVE:: loc_zspan TYPE(mumps_mat), SAVE :: femat !< Finite Element Method matrix for the full domain TYPE(mumps_mat), SAVE :: reduccedmat !< Finite Element Method matrix in the redduced web-spline sub-space !TYPE(mumps_mat), SAVE :: fematmpi !< Finite Element Method matrix prepared for mpi parallelism INTEGER :: nbmoments = 10 !< number of moments to be calculated and stored INTEGER(kind=omp_lock_kind), Allocatable:: mu_lock(:) !< Stores the lock for fields parallelism CONTAINS SUBROUTINE mag_init USE basic, ONLY: magnetfile, nr, nz USE bsplines USE mumps_bsplines USE mpihelper USE geometry ALLOCATE (Br((nr + 1)*(nz + 1)), Bz((nr + 1)*(nz + 1))) ALLOCATE (Athet((nr + 1)*(nz + 1))) ! Calculate magnetic field mirror components in grid points (Davidson analytical formula employed) ! or load it from magnetfile if present CALL magnet(magnetfile) end subroutine mag_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for solving Poisson and computes the magnetic field on the grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_init USE basic, ONLY: pot, nlperiod, nrank, rhs, volume, rgrid USE bsplines USE geometry USE mumps_bsplines USE mpihelper INTEGER :: nrz(2), i, d2, k1, n1 ! Auxiliary vectors ALLOCATE(vec1((nz+1)*(nr+1)),vec2((nr+1)*(nz+1))) DO i=0,nr vec1(i*(nz+1)+1:(i+1)*(nz+1))=zgrid!(0:nz) vec2(i*(nz+1)+1:(i+1)*(nz+1))=rgrid(i) END DO ! Set up 2d spline splrz used in the FEM CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz, nlppform=nlppform, period=nlperiod) ! Set up 2d spline splrz_ext used in the FEM to calculate the external electric field and potential CALL set_spline(femorder, ngauss, zgrid, rgrid, splrz_ext, nlppform=nlppform, period=nlperiod) !Allocate the work buffer to calculate the ppform d2 = splrz%sp2%dim k1 = splrz%sp1%order n1 = splrz%sp1%nints ALLOCATE(ppformwork(d2,k1,n1)) ! Calculate dimension of splines nrz(1) = nz nrz(2) = nr CALL get_dim(splrz, nrank, nrz, femorder) ! Allocate necessary variables ALLOCATE (matcoef(nrank(1), nrank(2))) ALLOCATE (pot((nr + 1)*(nz + 1))) ALLOCATE (potxt((nr + 1)*(nz + 1))) ALLOCATE (Erxt((nr + 1)*(nz + 1))) ALLOCATE (Ezxt((nr + 1)*(nz + 1))) ALLOCATE (rhs(nrank(1)*nrank(2))) ALLOCATE (gradgtilde(nrank(1)*nrank(2))) gradgtilde = 0 ALLOCATE (phi_spline(nrank(1)*nrank(2))) ALLOCATE (volume(nrank(1)*nrank(2))) volume = 0 ALLOCATE (Er((nr + 1)*(nz + 1)), Ez((nr + 1)*(nz + 1))) ALLOCATE (mu_lock(nrank(1)*nrank(2))) do i = 1, nrank(1)*nrank(2) call omp_init_lock(mu_lock(i)) end do end SUBROUTINE fields_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the geometry definition and read it from the standard input !> Precomputes the LHS matrix to solve Poisson abd the RHS effect of the dirichlet boundaries ! !--------------------------------------------------------------------------- SUBROUTINE fields_start USE geometry USE basic, ONLY: nrank implicit none INTEGER:: i,j, ierr DOUBLE PRECISION:: val ! set up the geometry module for setting up non-conforming boundary conditions call timera(0, "geom_init") call geom_init(splrz, vec1, vec2) call timera(1, "geom_init") ! Initialisation of FEM matrix CALL init(nrank(1)*nrank(2), 2, femat) ! Calculate and factorise FEM matrix (depends only on mesh) CALL fematrix(femat) If (walltype .lt. 0) then allocate (fverif(nrank(1)*nrank(2))) fverif = 0 end if ! Compute the volume of the splines and gtilde for solving E using web-splines CALL comp_volume !$OMP PARALLEL Call comp_gradgtilde !$OMP END PARALLEL if (nlweb) then ! Calculate reduced matrix for use of web splines call timera(0, "reduce femat") call Reducematrix(femat, reduccedmat) call timera(1, "reduce femat") call factor(reduccedmat) else call factor(femat) end if !WRITE(*,*) "Copy and to_mat worked" !CALL MPI_abort(MPI_COMM_WORLD,-1,ierr) !$OMP PARALLEL call vacuum_field !$OMP END PARALLEL END SUBROUTINE fields_start !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Recomputes the vacuum electric field ! !--------------------------------------------------------------------------- subroutine vacuum_field Use geometry USE basic, ONLY: pot, rhs implicit none + INTEGER:: i ! Computes the externally imposed electric field - !$OMP SINGLE - rhs = -gradgtilde - if (walltype .lt. 0) rhs = rhs + fverif - !$OMP END SINGLE + !$OMP DO SIMD + do i=1,nrank(1)*nrank(2) + rhs(i)=-gradgtilde(i) + !rhs = -gradgtilde + if (walltype .lt. 0) rhs(i) = rhs(i) + fverif(i) + end do + !$OMP END DO SIMD !$OMP BARRIER call poisson(splrz_ext) !$OMP BARRIER !$OMP SINGLE potxt = pot erxt = Er Ezxt = Ez !$OMP END SINGLE NOWAIT end subroutine !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Set-up the necessary variables for the communication of moments and rhs grid ! !--------------------------------------------------------------------------- SUBROUTINE fields_comm_init(Zbounds) USE basic, ONLY: nrank USE mpihelper INTEGER:: Zbounds(0:) loc_zspan = Zbounds(mpirank + 1) - Zbounds(mpirank) + femorder(1) if (allocated(loc_moments)) deallocate (loc_moments) ALLOCATE (loc_moments(nbmoments, loc_zspan*nrank(2))) if (allocated(loc_rhs)) deallocate (loc_rhs) ALLOCATE (loc_rhs(loc_zspan*nrank(2))) IF (mpisize .gt. 1) THEN CALL init_overlaps(nrank, femorder, Zbounds(mpirank), Zbounds(mpirank + 1), nbmoments) END IF END SUBROUTINE fields_comm_init !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Construct the right hand side vector used in the FEM Poisson solver ! !> @param[in] plist list of the particles type storing the desired specie parameters ! !--------------------------------------------------------------------------- SUBROUTINE rhscon(plist) USE bsplines use mpi USE basic, ONLY: rhs USE beam, ONLY: particles USE mpihelper Use geometry Use omp_lib type(particles), INTENT(INOUT):: plist(:) INTEGER:: i IF (nlphis) then ! We calculate the self-consistent field !$OMP SINGLE loc_rhs = 0 ! Reset the moments matrix !$OMP END SINGLE NOWAIT ! Assemble rhs for each specie Do i = 1, size(plist, 1) if (plist(i)%is_field) CALL deposit_charge(plist(i), loc_rhs) END Do !$OMP BARRIER ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL rhs_gather(rhs) ELSE !$OMP SINGLE rhs = loc_rhs !$OMP END SINGLE END IF ELSE ! We only consider the externally imposed field !$OMP SINGLE rhs = 0 !$OMP END SINGLE END IF !$OMP SINGLE IF (mpirank .eq. 0) THEN rhs = rhs - gradgtilde if (walltype .lt. 0) rhs = rhs + fverif END IF !$OMP END SINGLE END SUBROUTINE rhscon !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculate the 0th 1st and 2nd order moments of the particle p and stores it in moment ! !> @param[in] p the particles type storing the desired specie parameters !> @param[out] moment the 2d array storing the calculated moments ! !--------------------------------------------------------------------------- SUBROUTINE momentsdiag(p) USE bsplines use mpi USE beam, ONLY: particles USE mpihelper Use geometry type(particles), INTENT(INOUT):: p !REAL(kind=db), INTENT(INOUT):: moment(:, :) !$OMP SINGLE loc_moments = 0 ! Reset the moments matrix ! Assemble rhs !$OMP END SINGLE IF (p%Nploc .ne. 0) THEN CALL deposit_moments(p, loc_moments) END IF !$OMP SINGLE if(.not. allocated(p%moments))THEN if(mpirank.eq.0)THEN Allocate(p%moments(nbmoments,nrank(1)*nrank(2))) else Allocate(p%moments(0,0)) end if end if !$OMP END SINGLE ! If we are using MPI parallelism, reduce the rhs on the root process IF (mpisize .gt. 1) THEN CALL moments_gather(p%moments) ELSE !$OMP SINGLE p%moments = loc_moments !$OMP END SINGLE NOWAIT END IF END SUBROUTINE momentsdiag !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles moments (n,v,v^2) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_moments(p, p_loc_moments) USE bsplines use mpi USE basic, ONLY: Zbounds USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: p_loc_moments REAL(kind=db), DIMENSION(:, :), Allocatable:: omp_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db) :: vr, vthet, vz, coeff REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 64) ! Particle bunch size used when calling basfun ! Assemble rhs IF (p%Nploc .gt. 0) THEN !!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(zleft,rleft,jw,it,iend,irow,jcol,mu,k,vr,vz,vthet,coeff,fun,fun2) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node allocate(omp_loc_moments(size(p_loc_moments,1),size(p_loc_moments,2))) omp_loc_moments=0 !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle !CALL locintv(splrz%sp2, p%R(i:iend), rleft(1:k)) !CALL locintv(splrz%sp1, p%Z(i:iend), zleft(1:k)) rleft(1:k) = p%rindex(i:iend) zleft(1:k) = p%zindex(i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(3,i:iend), splrz%sp1, fun(:, :, 1:k), zleft(1:k) + 1) CALL basfun(p%pos(1,i:iend), splrz%sp2, fun2(:, :, 1:k), rleft(1:k) + 1) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) DO it = 1, (femorder(1) + 1) irow = zleft(k) + it - Zbounds(mpirank) jcol = rleft(k) + jw mu = irow + (jcol - 1)*(loc_zspan) coeff = p%weight*fun(it, 0, k)*fun2(jw, 0, k) ! Add contribution of particle nbunch to rhs grid point mu vr = 0.5*(p%U(1,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(1,i + k - 1)/p%Gammaold(i + k - 1)) vz = 0.5*(p%U(3,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(3,i + k - 1)/p%Gammaold(i + k - 1)) vthet = 0.5*(p%U(2,i + k - 1)/p%Gamma(i + k - 1) + p%Uold(2,i + k - 1)/p%Gammaold(i + k - 1)) !call omp_set_lock(mu_lock(mu)) !!$OMP ATOMIC UPDATE omp_loc_moments(1, mu) = omp_loc_moments(1, mu) + coeff !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(2, mu) = omp_loc_moments(2, mu) + coeff*vr !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(3, mu) = omp_loc_moments(3, mu) + coeff*vthet !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(4, mu) = omp_loc_moments(4, mu) + coeff*vz !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(5, mu) = omp_loc_moments(5, mu) + coeff*vr*vr !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(6, mu) = omp_loc_moments(6, mu) + coeff*vr*vthet !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(7, mu) = omp_loc_moments(7, mu) + coeff*vr*vz !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(8, mu) = omp_loc_moments(8, mu) + coeff*vthet*vthet !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(9, mu) = omp_loc_moments(9, mu) + coeff*vthet*vz !!$OMP END ATOMIC !!$OMP ATOMIC UPDATE omp_loc_moments(10, mu) = omp_loc_moments(10, mu) + coeff*vz*vz !!$OMP END ATOMIC !call omp_unset_lock(mu_lock(mu)) END DO END DO END DO END DO !!$OMP END PARALLEL DO !$OMP END DO NOWAIT Do i=1,size(p_loc_moments,2) call omp_set_lock(mu_lock(i)) p_loc_moments(:,i)=p_loc_moments(:,i)+omp_loc_moments(:,i) call omp_unset_lock(mu_lock(i)) end do !!$OMP END CRITICAL(loc_moments_reduce) DEALLOCATE (fun, fun2, zleft, rleft) END IF END subroutine deposit_moments !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Deposit the particles charges (q) from p on the grid ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] p_loc_moments local tensor used to store the moments of the given specie !--------------------------------------------------------------------------- SUBROUTINE deposit_charge(p, p_loc_moments) USE bsplines use mpi USE constants USE basic, ONLY: Zbounds, rnorm, phinorm USE beam, ONLY: particles USE mpihelper USE geometry USE omp_lib TYPE(particles), INTENT(IN):: p REAL(kind=db), DIMENSION(:), INTENT(INOUT):: p_loc_moments REAL(kind=db), DIMENSION(:), allocatable:: omp_loc_moments INTEGER ::irow, jcol, it, jw, mu, i, k, iend, nbunch INTEGER, DIMENSION(:), ALLOCATABLE::zleft, rleft REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) INTEGER:: num_threads, curr_thread real(kind=db):: contrib, chargecoeff num_threads = omp_get_max_threads() nbunch = p%Nploc/num_threads ! Particle bunch size used when calling basfun nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = min(nbunch, 16) ! Particle bunch size used when calling basfun chargecoeff = p%weight*p%q/(2*pi*eps_0*phinorm*rnorm) ! Normalized charge density simulated by each macro particle ! Assemble rhs IF (p%Nploc .ne. 0) THEN !!!$OMP PARALLEL DEFAULT(SHARED), PRIVATE(i,zleft, rleft, jw, it, iend, irow, jcol, mu, k, fun, fun2, contrib) ALLOCATE (zleft(nbunch), rleft(nbunch)) ALLOCATE (fun(1:femorder(1) + 1, 0:0, nbunch), fun2(1:femorder(2) + 1, 0:0, nbunch)) ! Arrays keeping values of b-splines at gauss node allocate(omp_loc_moments(size(p_loc_moments))) omp_loc_moments=0 zleft=0 rleft=0 curr_thread=omp_get_thread_num() !$OMP DO DO i = 1, p%Nploc, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, p%Nploc) k = iend - i + 1 ! Localize the particle rleft(1:k) = p%rindex(i:iend) zleft(1:k) = p%zindex(i:iend) ! Compute the value of the splines at the particles positions CALL basfun(p%pos(3,i:iend), splrz%sp1, fun, zleft(1:k) + 1) CALL basfun(p%pos(1,i:iend), splrz%sp2, fun2, rleft(1:k) + 1) !CALL geom_weight(p%Z(i:iend),p%R(i:iend),wgeom) DO k = 1, (iend - i + 1) DO jw = 1, (femorder(2) + 1) DO it = 1, (femorder(1) + 1) irow = zleft(k) + it - Zbounds(mpirank) jcol = rleft(k) + jw mu = irow + (jcol - 1)*(loc_zspan) ! Add contribution of particle k to rhs grid point mu contrib = fun(it, 0, k)*fun2(jw, 0, k)*p%geomweight(0,i + k - 1)*chargecoeff omp_loc_moments(mu) = omp_loc_moments(mu) + contrib END DO END DO END DO END DO !$OMP END DO NOWAIT DEALLOCATE (fun, fun2, zleft, rleft) !!!$OMP END PARALLEL !!$OMP CRITICAL(loc_charge_reduce) Do i=1,size(p_loc_moments) !$OMP ATOMIC p_loc_moments(i)=p_loc_moments(i)+omp_loc_moments(i) !$OMP END ATOMIC end do !!$OMP END CRITICAL(loc_charge_reduce) END IF END subroutine deposit_charge !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> and reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE rhs_gather(rhs) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:), INTENT(INOUT):: rhs INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) INTEGER:: overlap_type INTEGER:: rcvoverlap_type displs = Zbounds(0:mpisize - 1) counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1) counts(mpisize) = counts(mpisize) + femorder(1) !$OMP MASTER !WRITE(*,*) mpirank, "wE communicate overlap rhs" CALL rhsoverlapcomm(mpirank, leftproc, rightproc, loc_rhs, nrank, femorder, loc_zspan - femorder(1)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !$OMP DO SIMD DO j = 1, femorder(1) DO i = 1, nrank(2) loc_rhs((i - 1)*loc_zspan + j) = loc_rhs((i - 1)*loc_zspan + j)& & + rhsoverlap_buffer(nrank(2)*(j - 1) + i) END DO END DO !$OMP END DO SIMD END IF ! Set communication vector type overlap_type = rhsoverlap_type rcvoverlap_type = rcvrhsoverlap_type !$OMP MASTER IF (mpirank .eq. 0) THEN rhs = 0 END IF CALL MPI_GATHERV(loc_rhs, counts(mpirank + 1), rhsoverlap_type, & & rhs, counts, displs, rcvrhsoverlap_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE rhs_gather !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Do the communication of the local moment matrices between mpi workers for the overlap grid points !> and reduce the result on the host ! !--------------------------------------------------------------------------- SUBROUTINE moments_gather(moment) USE mpihelper USE Basic, ONLY: Zbounds, mpirank, leftproc, rightproc REAL(kind=db), DIMENSION(:, :), INTENT(INOUT):: moment INTEGER:: ierr, i, j INTEGER:: displs(mpisize), counts(mpisize) displs = Zbounds(0:mpisize - 1) counts = Zbounds(1:mpisize) - Zbounds(0:mpisize - 1) counts(mpisize) = counts(mpisize) + femorder(1) !$OMP MASTER CALL momentsoverlapcomm(mpirank, leftproc, rightproc, loc_moments, nrank, femorder, loc_zspan - femorder(1)) !$OMP END MASTER !$OMP BARRIER IF (mpirank .gt. 0) THEN !!$OMP PARALLEL DO SIMD DEFAULT(SHARED) private(i) !$OMP DO SIMD DO j = 1, femorder(1) DO i = 1, nrank(2) loc_moments(1:nbmoments, (i - 1)*loc_zspan + j) = loc_moments(1:nbmoments, (i - 1)*loc_zspan + j)& & + momentsoverlap_buffer(nbmoments*(nrank(2)*(j - 1) + i - 1) + 1:nbmoments*(nrank(2)*(j - 1) + i)) END DO END DO !$OMP END DO SIMD END IF !$OMP MASTER ! Set communication vector type IF (mpirank .eq. 0) THEN moment = 0 END IF CALL MPI_GATHERV(loc_moments, counts(mpirank + 1), momentsoverlap_type, & & moment, counts, displs, rcvmomentsoverlap_type, 0, MPI_COMM_WORLD, ierr) !$OMP END MASTER !$OMP BARRIER END SUBROUTINE moments_gather !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Solves Poisson equation using FEM. Distributes the result on all MPI workers and interpolate the electric forces !> for each particle. ! !--------------------------------------------------------------------------- SUBROUTINE poisson(splinevar) USE basic, ONLY: rhs, nrank, pot, nlend USE bsplines, ONLY: spline2d, gridval USE mumps_bsplines, ONLY: bsolve, vmx USE futils Use geometry type(spline2d):: splinevar INTEGER:: ierr real(kind=db), allocatable::reducedrhs(:) real(kind=db), allocatable:: reducedsol(:), tempcol(:) !$OMP MASTER if (nlweb) then ! we use the web-spline reduction for stability allocate (reducedrhs(nrank(1)*nrank(2))) allocate (reducedsol(nbreducedspline)) allocate (tempcol(nrank(1)*nrank(2))) if(mpirank.eq.0) then ! Only the root process solves Poisson reducedrhs = vmx(etilde, rhs) Call bsolve(reduccedmat, reducedrhs(1:nbreducedspline), reducedsol) end if CALL MPI_Bcast(reducedsol, nbreducedspline, db_type, 0, MPI_COMM_WORLD, ierr) tempcol = 0 tempcol(1:nbreducedspline) = reducedsol !phi_spline = 0 phi_spline = vmx(etildet, tempcol) else if(mpirank.eq.0) then CALL bsolve(femat, rhs, phi_spline) end if CALL MPI_Bcast(phi_spline, nrank(1)*nrank(2), db_type, 0, MPI_COMM_WORLD, ierr) end if matcoef = reshape(phi_spline, (/nrank(1), nrank(2)/)) !$OMP END MASTER !$OMP BARRIER ! update the ppform coefficients CALL updt_ppform2d(splinevar, matcoef) !$OMP BARRIER !$OMP SINGLE IF (mpirank .eq. 0 .and. (modulo(step, it2d) .eq. 0 .or. nlend)) THEN ! On the root process, compute the electric field for diagnostic purposes CALL gridval(splinevar, vec1, vec2, pot, (/0, 0/)) CALL gridval(splinevar, vec1, vec2, Ez, (/1, 0/)) CALL gridval(splinevar, vec1, vec2, Er, (/0, 1/)) Ez = -pot*gridwdir(1,:) - Ez*gridwdir(0,:) - gtilde(1,:) Er = -pot*gridwdir(2,:) - Er*gridwdir(0,:) - gtilde(2,:) pot = pot*gridwdir(0,:) + gtilde(0,:) END IF !$OMP END SINGLE NOWAIT END SUBROUTINE poisson !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the electric fields and potential at the particles position for particles !> between positions nstart and nend in the list ! !> @param[in] p the particles type storing the desired specie parameters !> @param[in] nstart starting index for the particle list !> @param[in] nend ending index for the particle list !--------------------------------------------------------------------------- SUBROUTINE EFieldscompatparts(p, nstart, nend) Use beam, ONLY: particles Use geometry Use splinebound TYPE(particles), INTENT(INOUT):: p INTEGER, OPTIONAL::nstart, nend INTEGER:: i, iend, nst, nnd INTEGER:: nbunch INTEGER:: num_threads Real(kind=db), ALLOCATABLE:: erext(:), ezext(:), gtildeloc(:, :) if (.not. present(nstart)) nst = 1 if (.not. present(nend)) nnd = p%Nploc !num_threads = omp_get_max_threads() !nbunch = (nnd - nst + 1)/num_threads ! Particle bunch size used when calling basfun !nbunch = max(nbunch, 1) ! Particle bunch size used when calling basfun nbunch = 64 ! Particle bunch size used when calling basfun Allocate (erext(nbunch), ezext(nbunch), gtildeloc(0:2,0:nbunch - 1)) ! Evaluate the electric potential and field at the particles position !$OMP DO DO i = nst, nnd, nbunch ! Avoid segmentation fault by accessing non relevant data iend = min(i + nbunch - 1, nnd) CALL speval(splrz, p%pos(3,i:iend), p%pos(1,i:iend),p%Zindex(i:iend),p%Rindex(i:iend), p%pot(i:iend), p%E(2,i:iend), p%E(1,i:iend)) CALL speval(splrz_ext, p%pos(3,i:iend), p%pos(1,i:iend),p%Zindex(i:iend),p%Rindex(i:iend), p%potxt(i:iend)) Call total_gtilde(p%pos(3,i:iend), p%pos(1,i:iend), gtildeloc(:,0:iend - i),p%geomweight(:,i:iend)) p%E(2,i:iend) = -p%E(2,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(1,i:iend) - gtildeloc(1,0:iend - i) p%E(1,i:iend) = -p%E(1,i:iend)*p%geomweight(0,i:iend) - p%pot(i:iend)*p%geomweight(2,i:iend) - gtildeloc(2,0:iend - i) p%pot(i:iend) = p%geomweight(0,i:iend)*p%pot(i:iend) + gtildeloc(0,0:iend - i) p%potxt(i:iend) = p%geomweight(0,i:iend)*p%potxt(i:iend) + gtildeloc(0,0:iend - i) END DO !$OMP END DO NOWAIT END SUBROUTINE EFieldscompatparts !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Constucts the FEM matrix using bsplines initialized in fields_init !--------------------------------------------------------------------------- SUBROUTINE fematrix(mat) USE bsplines USE geometry USE omp_lib USE sparse type(mumps_mat):: mat REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) - REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :) + REAL(kind=db), ALLOCATABLE :: fun(:, :, :), fun2(:, :, :) REAL(kind=db) :: contrib INTEGER, ALLOCATABLE :: idert(:, :), iderw(:, :), iderg(:, :) + integer,allocatable:: iid(:),jid(:) INTEGER :: i, j, jt, iw, irow, jcol, mu, igauss, iterm, irow2, jcol2, mu2, kterms, gausssize kterms=8 - ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node + + If (allocated(fun)) deallocate (fun) + If (allocated(fun2)) deallocate (fun2) + ALLOCATE (fun(1:femorder(1) + 1, 0:1,3*ngauss(1)*ngauss(2)), fun2(1:femorder(2) + 1, 0:1,3*ngauss(1)*ngauss(2))) + If (allocated(wgeom)) deallocate (wgeom) + ALLOCATE (wgeom(0:2,3*ngauss(1)*ngauss(2)))!Arrays keeping values of b-splines at gauss node + !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines ALLOCATE (idert(kterms, 2), iderw(kterms, 2), coefs(kterms), iderg(kterms, 2)) + ALLOCATE (iid(3*ngauss(1)*ngauss(2)), jid(3*ngauss(1)*ngauss(2))) !Pointers on the order of derivatives call timera(0, "fematrix") ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO ! Assemble FEM matrix -!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, iterm,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, iderw, idert, iderg, coefs, fun, fun2) +!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, iterm,jt,irow,jcol, mu, iw, irow2,jcol2, mu2, contrib, iderw, idert, iderg, coefs, fun, fun2,iid,jid) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position !! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) - if (gausssize .gt. 0) then - If (allocated(wgeom)) deallocate (wgeom) - ALLOCATE (wgeom(0:2,gausssize)) - CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom) + + iid=i + jid=j + + if (gausssize .gt. 1) then + !If (allocated(wgeom)) deallocate (wgeom) + !ALLOCATE (wgeom(0:2,gausssize)) + CALL geom_weight(xgauss(1:gausssize, 1), xgauss(1:gausssize, 2), wgeom(:,1:gausssize)) + CALL basfun(xgauss(1:gausssize, 1), splrz%sp1, fun(:,:,1:gausssize), iid(1:gausssize)) + CALL basfun(xgauss(1:gausssize, 2), splrz%sp2, fun2(:,:,1:gausssize), jid(1:gausssize)) End if - DO igauss = 1, gausssize ! Loop on gaussian weights and positions - CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i) - CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j) - CALL coefeq(xgauss(igauss, :), idert, iderw, iderg, coefs, kterms) - DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) - irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 - mu = irow + (jcol - 1)*nrank(1) - call omp_set_lock(mu_lock(mu)) - DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) - irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 - mu2 = irow2 + (jcol2 - 1)*nrank(1) - contrib=0.0_db + DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) + irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 + mu = irow + (jcol - 1)*nrank(1) + + DO iw = 1, (1 + femorder(1))*(femorder(2) + 1) + irow2 = i + f(iw, 1) - 1; jcol2 = j + f(iw, 2) - 1 + mu2 = irow2 + (jcol2 - 1)*nrank(1) + contrib=0.0_db + DO igauss = 1, gausssize ! Loop on gaussian weights and positions + CALL coefeq(xgauss(igauss, :), idert, iderw, iderg, coefs, kterms) DO iterm = 1, kterms ! Loop on the two integration dimensions contrib = contrib+wgeom(iderg(iterm, 1),igauss)*wgeom(iderg(iterm, 2),igauss)* & - & fun(f(jt, 1), idert(iterm, 1))*fun(f(iw, 1), idert(iterm, 2))* & - & fun2(f(jt, 2), iderw(iterm, 1))*fun2(f(iw, 2), iderw(iterm, 2))* & + & fun(f(jt, 1), idert(iterm, 1),igauss)*fun(f(iw, 1), idert(iterm, 2),igauss)* & + & fun2(f(jt, 2), iderw(iterm, 1),igauss)*fun2(f(iw, 2), iderw(iterm, 2),igauss)* & & wgauss(igauss)*coefs(iterm) END DO - CALL updt_sploc(mat%mat%row(mu), mu2, contrib) - END DO + end do + call omp_set_lock(mu_lock(mu)) + CALL updt_sploc(mat%mat%row(mu), mu2, contrib) call omp_unset_lock(mu_lock(mu)) END DO + END DO END DO END DO !$OMP End parallel do DEALLOCATE (f, aux) DEALLOCATE (idert, iderw, coefs, fun, fun2) call timera(1, "fematrix") END SUBROUTINE fematrix !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the volume of the splines cells needed to display the density in post-processing !--------------------------------------------------------------------------- SUBROUTINE comp_volume USE bsplines USE geometry USE basic, ONLY: Volume REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), gtildeintegr(:, :), ftestpt(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib call timera(0, "comp_volume") ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO volume = 0 if (walltype .lt. 0) fverif = 0 ! Assemble Volume matrix !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, gtildeintegr, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib) DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) If (allocated(wgeom)) deallocate (wgeom) if (gausssize .gt. 0) then ALLOCATE (wgeom(0:2,size(xgauss, 1))) CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom) End if If (allocated(gtildeintegr)) deallocate (gtildeintegr) ALLOCATE (gtildeintegr(0:2,size(xgauss, 1))) Call total_gtilde(xgauss(:, 1), xgauss(:, 2), gtildeintegr,wgeom) if (walltype .lt. 0) then If (allocated(ftestpt)) deallocate (ftestpt) ALLOCATE (ftestpt(0:0,size(xgauss, 1))) CALL ftest(xgauss(:, 1), xgauss(:, 2), ftestpt) end if DO igauss = 1, gausssize ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j) CALL coefeqext(xgauss(igauss, :), idt, idw, idg, idp, coefs) DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 2*pi*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)*wgauss(igauss)*xgauss(igauss, 2)!*wgeom(igauss,0) !$OMP ATOMIC UPDATE volume(mu) = volume(mu) + newcontrib !$OMP END ATOMIC if (walltype .lt. 0) THEN newcontrib = ftestpt(0,igauss)*fun(f(jt, 1), 0)*fun2(f(jt, 2), 0)& &*wgeom(0,igauss)*wgauss(igauss)*xgauss(igauss, 2) !$OMP ATOMIC UPDATE fverif(mu) = fverif(mu) + newcontrib !$OMP END ATOMIC end if END DO END DO END DO END DO !$OMP END PARALLEL DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) call timera(1, "comp_volume") END SUBROUTINE comp_volume !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the gradient of the gtilde function for the web-spline method needed to correctly apply the dirichlet boundary conditions !--------------------------------------------------------------------------- SUBROUTINE comp_gradgtilde USE bsplines USE geometry REAL(kind=db), ALLOCATABLE :: xgauss(:, :), wgauss(:), wgeom(:, :) INTEGER, ALLOCATABLE :: f(:, :), aux(:) REAL(kind=db), ALLOCATABLE :: coefs(:) REAL(kind=db), ALLOCATABLE :: fun(:, :), fun2(:, :), gtildeintegr(:, :), ftestpt(:, :) Integer, ALLOCATABLE, Dimension(:) :: idg, idt, idp, idw INTEGER :: i, j, jt, irow, jcol, mu, igauss, gausssize, iterm, nterms Real(kind=db)::newcontrib !call timera(0, "comp_gradgtilde") ALLOCATE (fun(1:femorder(1) + 1, 0:1), fun2(1:femorder(2) + 1, 0:1))!Arrays keeping values of b-splines at gauss node !ALLOCATE(xgauss(ngauss(1)*ngauss(2),2), wgauss(ngauss(1)*ngauss(2)),zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2))) !Gaussian nodes and weights arrays ALLOCATE (f((femorder(1) + 1)*(femorder(2) + 1), 2), aux(femorder(1) + 1)) !Auxiliary arrays ordering bsplines nterms = 4 Allocate (idg(nterms), idt(nterms), idw(nterms), idp(nterms), coefs(nterms)) ! Constuction of auxiliary array ordering bsplines in given interval DO i = 1, (femorder(1) + 1) aux(i) = i END DO DO i = 1, (femorder(2) + 1) f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 1) = aux f((i - 1)*(femorder(1) + 1) + 1:i*(femorder(1) + 1), 2) = i END DO !$OMP SINGLE gradgtilde = 0 !$OMP END SINGLE ! Assemble gradgtilde matrix !! $OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(j,i,xgauss,wgauss,gausssize,wgeom, igauss, gtildeintegr, ftestpt, iterm,jt,irow,jcol, mu, idw, idt, idg, idp, coefs, fun, fun2, newcontrib) !$OMP DO DO j = 1, nr ! Loop on r position DO i = 1, nz ! Loop on z position ! Computation of gauss weight and position in r and z direction for gaussian integration Call calc_gauss(splrz, ngauss, i, j, xgauss, wgauss, gausssize) If (allocated(wgeom)) deallocate (wgeom) if (gausssize .gt. 0) then ALLOCATE (wgeom(0:2,size(xgauss, 1))) CALL geom_weight(xgauss(:, 1), xgauss(:, 2), wgeom) End if If (allocated(gtildeintegr)) deallocate (gtildeintegr) ALLOCATE (gtildeintegr(0:2,size(xgauss, 1))) Call total_gtilde(xgauss(:, 1), xgauss(:, 2), gtildeintegr,wgeom) if (walltype .lt. 0) then If (allocated(ftestpt)) deallocate (ftestpt) ALLOCATE (ftestpt(0:0,size(xgauss, 1))) CALL ftest(xgauss(:, 1), xgauss(:, 2), ftestpt) end if DO igauss = 1, gausssize ! Loop on gaussian weights and positions CALL basfun(xgauss(igauss, 1), splrz%sp1, fun, i) CALL basfun(xgauss(igauss, 2), splrz%sp2, fun2, j) CALL coefeqext(xgauss(igauss, :), idt, idw, idg, idp, coefs) DO jt = 1, (1 + femorder(1))*(femorder(2) + 1) irow = i + f(jt, 1) - 1; jcol = j + f(jt, 2) - 1 mu = irow + (jcol - 1)*nrank(1) newcontrib = 0 Do iterm = 1, nterms newcontrib = newcontrib + wgeom( idg(iterm),igauss)*gtildeintegr( idp(iterm),igauss)* & & fun(f(jt, 1), idt(iterm))*fun2(f(jt, 2), idw(iterm))* & & wgauss(igauss)*coefs(iterm) End do !$OMP ATOMIC UPDATE gradgtilde(mu) = gradgtilde(mu) + newcontrib !$OMP END ATOMIC END DO END DO END DO END DO !!! $OMP END PARALLEL DO !$OMP END DO !DEALLOCATE(xgauss, wgauss,zg,rg, wzg, wrg) DEALLOCATE (f, aux) DEALLOCATE (fun, fun2) !call timera(1, "comp_gradgtilde") END SUBROUTINE comp_gradgtilde !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Imposes the dirichlet boundary conditions on the FEM matrix for the case where we use regular splines ( not web-splines). !--------------------------------------------------------------------------- SUBROUTINE fe_dirichlet REAL(kind=db), ALLOCATABLE :: arr(:) INTEGER :: i ALLOCATE (arr(nrank(1)*nrank(2))) DO i = 1, nrank(1) IF (rgrid(0) .ne. 0.0_db) THEN arr = 0; arr(i) = 1; CALL putrow(femat, i, arr) END IF arr = 0; arr(nrank(1)*nrank(2) + 1 - i) = 1; CALL putrow(femat, nrank(1)*nrank(2) + 1 - i, arr) END DO DEALLOCATE (arr) END SUBROUTINE fe_dirichlet !________________________________________________________________________________ SUBROUTINE coefeq(x, idt, idw, idg, c, kterms) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idt(:, :), idw(:, :), idg(:, :),kterms REAL(kind=db), INTENT(out) :: c(:) kterms=8 c = x(2) idt(1, 1) = 0 idt(1, 2) = 0 idw(1, 1) = 0 idw(1, 2) = 0 idg(1, 1) = 1 idg(1, 2) = 1 idt(2, 1) = 0 idt(2, 2) = 1 idw(2, 1) = 0 idw(2, 2) = 0 idg(2, 1) = 1 idg(2, 2) = 0 idt(3, 1) = 1 idt(3, 2) = 0 idw(3, 1) = 0 idw(3, 2) = 0 idg(3, 1) = 0 idg(3, 2) = 1 idt(4, 1) = 1 idt(4, 2) = 1 idw(4, 1) = 0 idw(4, 2) = 0 idg(4, 1) = 0 idg(4, 2) = 0 idt(5, 1) = 0 idt(5, 2) = 0 idw(5, 1) = 0 idw(5, 2) = 0 idg(5, 1) = 2 idg(5, 2) = 2 idt(6, 1) = 0 idt(6, 2) = 0 idw(6, 1) = 0 idw(6, 2) = 1 idg(6, 1) = 2 idg(6, 2) = 0 idt(7, 1) = 0 idt(7, 2) = 0 idw(7, 1) = 1 idw(7, 2) = 0 idg(7, 1) = 0 idg(7, 2) = 2 idt(8, 1) = 0 idt(8, 2) = 0 idw(8, 1) = 1 idw(8, 2) = 1 idg(8, 1) = 0 idg(8, 2) = 0 END SUBROUTINE coefeq SUBROUTINE coefeqext(x, idt, idw, idg, idp, c) REAL(kind=db), INTENT(in) :: x(:) INTEGER, INTENT(out) :: idp(:), idt(:), idw(:), idg(:) REAL(kind=db), INTENT(out) :: c(:) c(1) = x(2) idp(1) = 1 idg(1) = 1 idt(1) = 0 idw(1) = 0 c(2) = x(2) idp(2) = 1 idg(2) = 0 idt(2) = 1 idw(2) = 0 c(3) = x(2) idp(3) = 2 idg(3) = 2 idt(3) = 0 idw(3) = 0 c(4) = x(2) idp(4) = 2 idg(4) = 0 idt(4) = 0 idw(4) = 1 END SUBROUTINE coefeqext !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Computes the magnetic field on the grid according to a magnetic mirror, !> or according to the linear interpolation of the values on the !> grid saved in h5 file stored at magfile. !> @param[in] magfile filname of .h5 file containing the definitions of A and B !--------------------------------------------------------------------------- SUBROUTINE magnet(magfile) USE basic, ONLY: B0, Rcurv, rgrid, zgrid, width, rnorm, nr, nz, bnorm USE constants, ONLY: Pi CHARACTER(LEN=*), INTENT(IN), OPTIONAL:: magfile REAL(kind=db) :: rg, zg, halfLz, MirrorRatio INTEGER :: i, rindex IF (len_trim(magfile) .lt. 1) THEN halfLz = (zgrid(nz) + zgrid(0))/2 MirrorRatio = (Rcurv - 1)/(Rcurv + 1) DO i = 1, (nr + 1)*(nz + 1) rindex = (i - 1)/(nz + 1) rg = rgrid(rindex) zg = zgrid(i - rindex*(nz + 1) - 1) - halfLz Br(i) = -B0*MirrorRatio*SIN(2*pi*zg/width*rnorm)*bessi1(2*pi*rg/width*rnorm)/bnorm Bz(i) = B0*(1 - MirrorRatio*COS(2*pi*zg/width*rnorm)*bessi0(2*pi*rg/width*rnorm))/bnorm Athet(i) = 0.5*B0*(rg*rnorm - width/pi*MirrorRatio*bessi1(2*pi*rg/width*rnorm)*COS(2*pi*zg/width*rnorm)) END DO ELSE CALL load_mag_from_h5(magfile) END IF END SUBROUTINE magnet !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Loads the magnetic field defined in the .h5 file at location magfile !> @param[in] magfile filname of .h5 file containing the definitions of A and B !--------------------------------------------------------------------------- SUBROUTINE load_mag_from_h5(magfile) USE basic, ONLY: B0, rnorm, bnorm, bscaling USE constants, ONLY: Pi USE futils USE bsplines CHARACTER(LEN=*), INTENT(IN):: magfile REAL(kind=db), ALLOCATABLE :: magr(:), magz(:) REAL(kind=db), ALLOCATABLE :: tempBr(:, :), tempBz(:, :), tempAthet(:, :) real(kind=db), allocatable:: c(:,:) type(spline2d):: Maginterpolation REAL(kind=db) :: maxB INTEGER :: magfid, dims(2) LOGICAL:: B_is_saved INTEGER :: magn(2), magrank CALL openf(trim(magfile), magfid, 'r', real_prec='d') CALL getdims(magfid, '/mag/Athet', magrank, magn) ALLOCATE (magr(magn(2)), magz(magn(1))) ALLOCATE (tempAthet(magn(1), magn(2)), tempBr(magn(1), magn(2)), tempBz(magn(1), magn(2))) ! Read r and z coordinates for the definition of A_\thet, and B CALL getarr(magfid, '/mag/r', magr) CALL getarr(magfid, '/mag/z', magz) CALL getarr(magfid, '/mag/Athet', tempAthet) IF (isdataset(magfid, '/mag/Br') .and. isdataset(magfid, '/mag/Bz')) THEN CALL getarr(magfid, '/mag/Br', tempBr) CALL getarr(magfid, '/mag/Bz', tempBz) IF(bscaling .gt. 0) then maxB=sqrt(maxval(tempBr**2+tempBz**2)) tempBr=tempBr/maxB*B0 tempBz=tempBz/maxB*B0 end if B_is_saved = .true. ELSE B_is_saved = .false. END IF magz=magz/rnorm magr=magr/rnorm CALL set_splcoef((/3,3/),magz,magr,Maginterpolation) call get_dim(Maginterpolation,dims) ! Interpolation of the magnetic potential vector allocate(c(dims(1),dims(2))) call get_splcoef(Maginterpolation,tempAthet, c) CALL gridval(Maginterpolation,vec1,vec2, Athet ,(/0,0/),c) if(B_is_saved == .true.)then ! Interpolation of the Axial magnetic field call get_splcoef(Maginterpolation,tempBz, c) CALL gridval(Maginterpolation,vec1,vec2, Bz ,(/0,0/),c) ! Interpolation of the radial magnetic field call get_splcoef(Maginterpolation,tempBr, c) CALL gridval(Maginterpolation,vec1,vec2, Br ,(/0,0/),c) else CALL gridval(Maginterpolation,vec1,vec2, Br,(/1,0/)) Br=-Br CALL gridval(Maginterpolation,vec1,vec2, Bz,(/0,1/)) Bz=Bz+Athet/vec2 end if if( bscaling .lt. 0 ) then maxB = maxval(sqrt(Bz**2 + Br**2)) Bz = Bz/maxB*B0 Br = Br/maxB*B0 end if ! We normalize Br = Br/bnorm Bz = Bz/bnorm CALL closef(magfid) deallocate(c) call destroy_SP(Maginterpolation) END SUBROUTINE load_mag_from_h5 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the zero order FUNCTION bessi0(x) REAL(kind=db) :: bessi0, x REAL(kind=db) :: ax REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9 DATA p1, p2, p3, p4, p5, p6, p7/1.0d0, 3.5156229d0, 3.0899424d0, 1.2067492d0, 0.2659732d0, 0.360768d-1, 0.45813d-2/ DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, 0.1328592d-1, 0.225319d-2, -0.157565d-2, 0.916281d-2, & & -0.2057706d-1, 0.2635537d-1, -0.1647633d-1, 0.392377d-2/ if (abs(x) .lt. 3.75) then y = (x/3.75)**2 bessi0 = p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7))))) else ax = abs(x) y = 3.75/ax bessi0 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))) end if return END FUNCTION bessi0 !________________________________________________________________________________ !Modified Bessel functions of the first kind of the first order FUNCTION bessi1(x) REAL(kind=db) :: bessi1, x REAL(kind=db) :: ax REAL(kind=db) p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9, y SAVE p1, p2, p3, p4, p5, p6, p7, q1, q2, q3, q4, q5, q6, q7, q8, q9 DATA p1, p2, p3, p4, p5, p6, p7/0.5d0, 0.87890594d0, 0.51498869d0, 0.15084934d0, 0.2658733d-1, 0.301532d-2, 0.32411d-3/ DATA q1, q2, q3, q4, q5, q6, q7, q8, q9/0.39894228d0, -0.3988024d-1, -0.362018d-2, 0.163801d-2, -0.1031555d-1, & & 0.2282967d-1, -0.2895312d-1, 0.1787654d-1, -0.420059d-2/ if (abs(x) .lt. 3.75D0) then y = (x/3.75D0)**2 bessi1 = x*(p1 + y*(p2 + y*(p3 + y*(p4 + y*(p5 + y*(p6 + y*p7)))))) else ax = abs(x) y = 3.75D0/ax bessi1 = (exp(ax)/sqrt(ax))*(q1 + y*(q2 + y*(q3 + y*(q4 + y*(q5 + y*(q6 + y*(q7 + y*(q8 + y*q9)))))))) if (x .lt. 0.) bessi1 = -bessi1 end if return END FUNCTION bessi1 !--------------------------------------------------------------------------- !> @author !> Patryk kaminski EPFL/SPC !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Free the memory used by the fields module !--------------------------------------------------------------------------- SUBROUTINE clean_fields Use bsplines USE basic, ONLY: rhs INTEGER:: i do i = 1, nrank(1)*nrank(2) call omp_destroy_lock(mu_lock(i)) end do DEALLOCATE (mu_lock) DEALLOCATE (matcoef) DEALLOCATE (pot) DEALLOCATE (rhs) DEALLOCATE (loc_rhs) DEALLOCATE (loc_moments) DEALLOCATE (phi_spline) DEALLOCATE (Br, Bz) DEALLOCATE (Er, Ez) DEALLOCATE (vec1, vec2) Call DESTROY_SP(splrz) Call DESTROY_SP(splrz_ext) END SUBROUTINE clean_fields SUBROUTINE updt_sploc(arow, j, val) ! ! Update element j of row arow or insert it in an increasing "index" ! USE sparse TYPE(sprow), TARGET :: arow INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(in) :: val ! TYPE(elt), TARGET :: pre_root TYPE(elt), POINTER :: t, p ! pre_root%next => arow%row0 ! pre_root is linked to the head of the list. t => pre_root DO WHILE (ASSOCIATED(t%next)) p => t%next IF (p%index .EQ. j) THEN p%val = p%val + val RETURN END IF IF (p%index .GT. j) EXIT t => t%next END DO ALLOCATE (p) p = elt(j, val, t%next) t%next => p ! arow%nnz = arow%nnz + 1 arow%row0 => pre_root%next ! In case the head is altered END SUBROUTINE updt_sploc SUBROUTINE updt_ppform2d(sp,c) use bsplines TYPE(spline2d), INTENT(inout) :: sp DOUBLE PRECISION, DIMENSION(:,:), INTENT(in) :: c !DOUBLE PRECISION, ALLOCATABLE :: work(:,:,:) INTEGER:: m,mm INTEGER :: d1, d2, k1, k2, n1, n2 d1 = sp%sp1%dim d2 = sp%sp2%dim k1 = sp%sp1%order k2 = sp%sp2%order n1 = sp%sp1%nints n2 = sp%sp2%nints !ALLOCATE(work(d2,k1,n1)) !$OMP DO DO m=1,SIZE(c,2) CALL topp0(sp%sp1, c(:,m), ppformwork(m,:,:)) END DO !$OMP END DO NOWAIT !$OMP SINGLE IF( ASSOCIATED(sp%ppform) ) DEALLOCATE(sp%ppform) ALLOCATE(sp%ppform(k1,n1,k2,n2)) !$OMP END SINGLE !$OMP DO DO mm=1,SIZE(ppformwork,3) DO m=1,SIZE(ppformwork,2) CALL topp0(sp%sp2, ppformwork(:,m,mm), sp%ppform(m,mm,:,:)) END DO END DO !$OMP END DO !DEALLOCATE(work) end subroutine updt_ppform2d !=========================================================================== SUBROUTINE topp0(sp, c, ppform) ! ! Compute PPFORM of a fuction defined by the spline SP ! and spline coefficients C(1:d) ! use bsplines TYPE(spline1d), INTENT(in) :: sp DOUBLE PRECISION, INTENT(in) :: c(:) DOUBLE PRECISION, INTENT(out) :: ppform(0:,:) INTEGER :: p, nints, i, j, k ! p = sp%order - 1 nints = sp%nints ! ppform = 0.0d0 DO i=1,nints ! on each knot interval DO j=1,p+1 ! all spline in interval i DO k=0,p ! k_th derivatives ppform(k,i) = ppform(k,i) + sp%val0(k,j,i)*c(j+i-1) END DO END DO END DO ! END SUBROUTINE topp0 !+ END MODULE fields diff --git a/src/geometry_mod.f90 b/src/geometry_mod.f90 index 3db46a4..5ea2498 100644 --- a/src/geometry_mod.f90 +++ b/src/geometry_mod.f90 @@ -1,1235 +1,1235 @@ !------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: geometry ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for handling geometries with non constant radius using b-splines interpolation !> This module defines ways to comupte the weight function needed for weighted extended b-splines and !> can load the definition of the geometry from input files !> This module is based on the theory by K. Hollig and book "Finite element methods with b-splines" !> SIAM Frontiers in applied mathematics !------------------------------------------------------------------------------ MODULE geometry USE constants USE bsplines USE mumps_bsplines USE splinebound use weighttypes IMPLICIT NONE type innerspline Integer, Allocatable:: k(:) ! Index in reduced set real(kind=db), Allocatable:: weight(:) ! geomtric weight at relevant cell end type innerspline type test_params !< parameters defining the manufactured test solutionwhen negative weighttypes is used real(kind=db):: z0 real(kind=db):: r0 real(kind=db):: Lz real(kind=db):: Lr end type Integer, save:: testkr=1, testkz=1 Logical, save:: nlweb=.true. !< use weighted extended b-splines and not only weighted b-splies Integer, save :: walltype = 0 !< type of geometric weight to use (see readgeom) Integer, save, Allocatable:: bsplinetype(:) !< Array containing the inner/outer type for each bspline Integer, save, Allocatable:: gridcelltype(:,:) !< Array containing the inner/outer type for each gridcell Integer, save, Allocatable:: linkedspline(:,:) !< Array containing the lowerleft linked spline in case of boundary spline Real(kind=db), save, allocatable:: gridwdir(:,:) !< Stores the Dirichlet geometric weight at the grid points Real(kind=db), save, allocatable:: gridwdom(:,:) !< Stores the domain weight at the grid points Real(kind=db), save, allocatable:: gtilde(:,:) ! Stores the extension to the domain of the boundary conditions type(mumps_mat):: etilde !> Matrix of extendend web splines definition 4.9 p48 of Hollig's book type(mumps_mat):: etildet ! Transpose of Matrix of extendend web splines integer,save :: nbreducedspline ! Number of splines in the reduced set type(test_params), save:: test_pars PROCEDURE(geom_eval), POINTER:: dirichlet_weight => NULL()!< Function evaluating the weight for Dirichelt boundary conditions PROCEDURE(geomtot_eval), POINTER:: domain_weight => NULL() !< function giving the limits of the simulation domain PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() !< Computes the parameter gtilde used to impose dirichlet boundary conditions with phi=uh+gtilde PUBLIC:: geom_weight, dom_weight ABSTRACT INTERFACE SUBROUTINE gtilde_eval(z,r,g,w) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) END SUBROUTINE SUBROUTINE geom_eval(z,r,w,wupper) USE constants Real(kind=db), INTENT(IN) :: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), OPTIONAL :: wupper(0:,:) END SUBROUTINE SUBROUTINE geomtot_eval(z,r,w,idwall) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) INTEGER, optional, INTENT(OUT):: idwall(:) END SUBROUTINE END INTERFACE INTERFACE geom_weight MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2 END INTERFACE geom_weight INTERFACE dom_weight MODULE PROCEDURE dom_weight0, dom_weight1, dom_weight2, dom_weight3 END INTERFACE dom_weight NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, z_a, z_b ,walltype, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight, testkr, testkz contains ! Read the input parameters from the standard input file SUBROUTINE read_geom(Fileid, rnorm, splrz, Potinn, Potout) use mpi Use bsplines use basic, ONLY: phinorm use weighttypes type(spline2d):: splrz Real(kind=db):: rnorm ! normalisation variable for distances Real(kind=db):: Potinn ! potential at inner electrode from basic Real(kind=db):: Potout ! potential at outer electrode from basic Integer:: Fileid, mpirank, ierr, istat character(len=1000) :: line CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) Rewind(Fileid) READ(Fileid, geomparams, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(Fileid) read(Fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, geomparams) !! Normalizations and initialization of geometric variables r_a=r_a/rnorm r_b=r_b/rnorm z_a=z_a/rnorm z_b=z_b/rnorm r_bLeft=r_bLeft/rnorm r_bRight=r_bRight/rnorm if(r_a .eq. 0 .and. r_b .eq.0) then !! in case no geom_params have been definedwe take the defaults from the grid limits r_a=splrz%sp2%knots(0) r_b=splrz%sp2%knots(splrz%sp2%nints) end if z_0=z_0/rnorm r_0=r_0/rnorm z_r=z_r/rnorm r_r=r_r/rnorm invr_r=1/r_r invr_z=1/z_r Phidown=Potinn Phiup=Potout SELECT CASE (abs(walltype)) CASE (2) ! coaxial cylinder and top ellipse with cylinder extensions total_gtilde=>gUpDown Dirichlet_weight=>geom_w2 domain_weight=>geom_rvaschtot CASE (3) ! Two ellipses with "parallel" tangents with cylinders total_gtilde=>gUpDown Dirichlet_weight=>geom_w3 domain_weight=>geom_rvaschtot CASE (4) ! Two ellipses with same radii with cylinders total_gtilde=>gUpDown Dirichlet_weight=>geom_w4 domain_weight=>geom_rvaschtot CASE (5) ! Two ellipses with same radii with cylinders and total_gtilde=>gUpDown Dirichlet_weight=>geom_w5 domain_weight=>geom_rvaschtot CASE (6) ! circular coaxial tilted ellipse right Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w6 domain_weight=>geom_rvaschtot CASE (7) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w7 domain_weight=>geom_rvaschtot CASE (8) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w8 domain_weight=>geom_rvaschtot CASE (9) ! Geometry defined as a spline curve total_gtilde=>gspline Dirichlet_weight=>geom_spline domain_weight=>geom_splinetot call read_splinebound(Fileid,the_domain, splrz, rnorm, phinorm) CASE (10) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w10 CASE (11) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w11 CASE (12) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w12 CASE DEFAULT ! Ellipse as in gt170 standard weight and straight coaxial configs total_gtilde=>gstd Dirichlet_weight=>geom_weightstd domain_weight=>geom_rvaschtot END SELECT ! If we are lauching a test case, we load the test_pars variable used for the manufactured solution if(walltype.lt.0) then test_pars%Lr=(splrz%sp2%knots(splrz%sp2%nints)-splrz%sp2%knots(0))/testkr test_pars%Lz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/testkz test_pars%r0=0.5*(splrz%sp2%knots(splrz%sp2%nints)+splrz%sp2%knots(0)) test_pars%z0=0.5*(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0)) total_gtilde=>gtest end if end subroutine read_geom ! Initialises the module and precomputes the matrix e used for extende splines ! Classify every grid cell (inner, outer, boundary) SUBROUTINE geom_init(spl2, vec1, vec2) type(spline2d):: spl2 real(kind=db):: vec1(:),vec2(:) Real(kind=db), Allocatable:: zgrid(:),rgrid(:) Integer:: nrank(2), nrz(2) Call get_dim(spl2, nrank, nrz) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) ! create a table of the classification for each cell Allocate(gridcelltype(nrz(1),nrz(2))) Call classifycell(zgrid, rgrid, gridcelltype) if (nlweb) then ! if we use extended splines, we need to build the e matrix linking inner and outer splines Allocate(bsplinetype(nrank(1)*nrank(2))) Call classifyspline(spl2, gridcelltype, bsplinetype) Call buildetilde(spl2,bsplinetype,gridcelltype) end if ! Precompute the domain and dirichlet weights at the grid/cell positions ALLOCATE(gridwdir(0:2,(nrz(1)+1)*(nrz(2)+1))) gridwdir=0 ALLOCATE(gridwdom(0:2,(nrz(1)+1)*(nrz(2)+1))) gridwdom=0 ALLOCATE(gtilde(0:2,(nrz(1)+1)*(nrz(2)+1))) gtilde=0 Call geom_weight (vec1, vec2, gridwdir) Call dom_weight (vec1, vec2, gridwdom) ! Precompute the gtilde at the grid cell position CALL total_gtilde(vec1, vec2, gtilde, gridwdir) end Subroutine geom_init ! Save this module run parameters to the h5 result file in the correct group Subroutine geom_diag(File_handle, str, rnorm) use mpi Use futils use weighttypes Integer:: File_handle Real(kind=db):: rnorm Character(len=*):: str CHARACTER(len=128):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0) THEN Write(grpname,'(a,a)') trim(str),"/geometry" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "r_a", r_a*RNORM) Call attach(File_handle, trim(grpname), "r_b", r_b*RNORM) Call attach(File_handle, trim(grpname), "z_a", z_a*RNORM) Call attach(File_handle, trim(grpname), "z_b", z_b*RNORM) Call attach(File_handle, trim(grpname), "z_0", z_0*RNORM) Call attach(File_handle, trim(grpname), "r_0", r_0*RNORM) Call attach(File_handle, trim(grpname), "r_r", r_r*RNORM) Call attach(File_handle, trim(grpname), "z_r", z_r*RNORM) Call attach(File_handle, trim(grpname), "L_r", test_pars%Lr*RNORM) Call attach(File_handle, trim(grpname), "L_z", test_pars%Lz*RNORM) Call attach(File_handle, trim(grpname), "interior", interior) Call attach(File_handle, trim(grpname), "above1", above1) Call attach(File_handle, trim(grpname), "above2", above2) Call attach(File_handle, trim(grpname), "walltype", walltype) Call putarr(File_handle, trim(grpname)//'/geomweight',transpose(gridwdom)) Call putarr(File_handle, trim(grpname)//'/dirichletweight',transpose(gridwdir)) Call putarr(File_handle, trim(grpname)//'/gtilde',transpose(gtilde)) Call putarr(File_handle, trim(grpname)//'/ctype',gridcelltype) Call putarr(File_handle, trim(grpname)//'/linked_s',linkedspline) Call putarr(File_handle, trim(grpname)//'/bsplinetype',bsplinetype) END IF End subroutine geom_diag ! Construct the e matrix used to link inner and outer b-splines Subroutine buildetilde(spl2, bsplinetype, celltype) USE mumps_bsplines Use basic, ONLY: rnorm, mpirank type(spline2d):: spl2 Integer:: bsplinetype(:) Integer:: celltype(1:,1:) Logical, Allocatable:: Ibsplinetype(:) Integer:: i,j,k, icellz, icellr, jcellz, jcellr, nrank(2), nrz(2), norder(2) real(kind=db), allocatable:: zgrid(:), rgrid(:) real(kind=db):: wgeomi, indexdistance, eij integer:: l,m, n, linkedi type(innerspline):: innersplinelist real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:) integer, allocatable:: igridmesh(:), jgridmesh(:) !if(allocated(etilde)) deallocate(etilde) nbreducedspline=count(bsplinetype .eq. 1) Call get_dim(spl2,nrank,nrz,norder) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) allocate(linkedspline(nrank(1),nrank(2))) allocate(innersplinelist%k(nrank(1)*nrank(2)),innersplinelist%weight(nrank(1)*nrank(2))) allocate(Ibsplinetype(nrank(1)*nrank(2))) allocate(rgridmesh(nrank(1)*nrank(2)), zgridmesh(nrank(1)*nrank(2))) allocate(igridmesh(nrank(1)*nrank(2)), jgridmesh(nrank(1)*nrank(2))) allocate(hmesh(nrank(1)*nrank(2))) allocate(d(size(bsplinetype))) Ibsplinetype=.False. ! Compute the center of the 2D splines for the Lagrange interpolation call calcsplinecenters(spl2%sp1,cz,hz) call calcsplinecenters(spl2%sp2,cr,hr) Do i=0,nrank(2)-1 zgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cz rgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cr(i+1) hmesh(i*nrank(1)+1:(i+1)*nrank(1))=hr(i+1)*hz igridmesh(i*nrank(1)+1:(i+1)*nrank(1))=(/ (j,j=0,nrank(1)-1)/) jgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=i End do ! allocate memory for etilde and its transpose call init(nrank(1)*nrank(2),nbreducedspline,etilde) call init(nrank(1)*nrank(2),nbreducedspline,etildet) k=1 ! Compute the terms eii for the inner b-splines Do i=1,nrank(1)*nrank(2) if(bsplinetype(i) .ne. 1) cycle ! span of this spline is completely outside D ! one cell of bspline i is completely in D icellz=mod(i-1,nrank(1))+1 icellr=(i-1)/(nrank(1))+1 outer: do l=max(1,icellz-norder(1)),min(nrz(1),icellz) do m=max(1,icellr-norder(2)),min(nrz(2),icellr) if (celltype(l,m).eq.1) then call geom_weight((zgrid(l)+zgrid(l+1))/2,(rgrid(m)+rgrid(m+1))/2,wgeomi) EXIT outer end if end do end do outer call putele(etilde,k,i,1/wgeomi) call putele(etildet,i,k,1/wgeomi) innersplinelist%k(i)=k innersplinelist%weight(i)=1/wgeomi k=k+1 linkedspline(icellz,icellr)=i Ibsplinetype(i)=icellz+norder(1) .le. nrank(1) .and. icellr+norder(2) .le. nrank(2) if(.not. Ibsplinetype(i)) cycle do m=0,norder(2) do l=0,norder(1) ! Check if all positive splines in this spline domain are inner splines Ibsplinetype(i)=Ibsplinetype(i) .and. (bsplinetype(i+l+m*nrank(1)) .eq. 1) end do end do end do ! Compute the terms eij for the outer b-splines !!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(d,k,icellz,icellr,jcellr,jcellz,indexdistance,n,i,m,l,eij,linkedi) Do j=1,nrank(1)*nrank(2) ! find the closest b-spline fully in D for the web method if(bsplinetype(j) .ne. 0) cycle d=0 where(Ibsplinetype)! calculate distance between center of Interior splines and spline j d=(zgridmesh(j)-zgridmesh)**2+(rgridmesh(j)-rgridmesh)**2 end where k=minloc(d,1,MASK=Ibsplinetype) icellz=mod(k-1,nrank(1)) icellr=(k-1)/(nrank(1)) jcellz=mod(j-1,nrank(1)) jcellr=(j-1)/(nrank(1)) indexdistance=real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db) if(d(k) .gt. ((norder(1)+2)**2+(norder(2)+2)**2)*2 .and. mpirank .eq. 0)then Write(*,'(a)') 'Warning on system conditioning, the number of radial or axial points could be too low!' Write(*,'(a,1f6.2,a,2(1pe12.4))') 'Distance found: ', sqrt(indexdistance), ' at (z,r): ',zgridmesh(j)*rnorm, rgridmesh(j)*rnorm !stop end if ! Compute the Lagrange polynomia linking spline i and spline j linkedspline(jcellz+1,jcellr+1)=k do n=0,norder(2) do i=0,norder(1) eij=1 !eij=1 do m=0,norder(2) if(n.eq.m) cycle eij=eij*(rgridmesh(j)-rgridmesh(k+m*nrank(1)))/(rgridmesh(k+n*nrank(1))-rgridmesh(k+m*nrank(1))) !eij=eij*real((jcellr-icellr-m),db)/real((n-m),db) end do do l=0,norder(1) if( i .eq. l ) cycle eij=eij*(zgridmesh(j)-zgridmesh(k+l))/(zgridmesh(k+i)-zgridmesh(k+l)) !eij=eij*real((jcellz-icellz-l),db)/real((i-l),db) end do linkedi=innersplinelist%k(k+i+n*nrank(1)) ! equivalent to findloc, necessary for ifort 17 eij=eij*innersplinelist%weight(k+i+n*nrank(1)) ! add the polynomia to the etilde matrix call putele(etilde,linkedi,j,eij) call putele(etildet,j,linkedi,eij) end do end do end do !!$OMP End parallel do call to_mat(etilde) call to_mat(etildet) end subroutine ! Routine to compute the center of the 2d b-splines used in Lagrange interpolation Subroutine calcsplinecenters(spl,ctrs,heights) use bsplines type(spline1d):: spl real(kind=db), allocatable:: ctrs(:), heights(:) integer:: nrank, nx, order, i, left1, j, left2, left3 real(kind=db):: x1, x2, x3 real(kind=db), allocatable:: fun1(:,:), fun2(:,:), fun3(:,:) real(kind=db),allocatable:: init_heights(:) call get_dim(spl, nrank, nx, order) if (allocated(ctrs)) deallocate(ctrs) if (allocated(heights)) deallocate(heights) allocate(heights(nrank)) allocate(ctrs(nrank)) allocate(fun1(1:order+1,0:1), fun2(1:order+1,0:1), fun3(1:order+1,0:1)) allocate(init_heights(nrank)) init_heights=1.0 ctrs(:)=spl%knots(0:nrank-1) call gridval(spl, ctrs, heights, 0, init_heights) if (order .gt.1) then do i=2,nrank-1 x1=spl%knots(i-order)+1e5*Epsilon(x1) x2=spl%knots(i-1)-1e5*Epsilon(x2) left1=max(0,i-order) left2=min(nx-2,i-2) call basfun(x1,spl,fun1, left1+1) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold call basfun(x2,spl,fun2, left2+1) fun3=1 j=0 Do while( j .le. 300) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold x3=(x1+x2)/2 call locintv(spl, x3, left3) call basfun(x3,spl,fun3, left3+1) if( (x2-x1).lt.1e-13) exit if(abs(fun3(i-left3,1)).lt.1e-15) exit if(fun3(i-left3,1)*fun1(i-left1,1).le.0) then fun2=fun3 x2=x3 left2=left3 else fun1=fun3 x1=x3 left1=left3 end if j=j+1 End do ctrs(i)=x3 heights(i)=fun3(i-left3,0) end do end if end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine classifycell(zgrid,rgrid,gridctype) Real(kind=db):: zgrid(:),rgrid(:) Integer:: gridctype(:,:) Integer::i,j gridctype=-1 !$OMP parallel do private(i) Do j=1,size(rgrid,1)-1 Do i=1,size(zgrid,1)-1 ! Determines the type inner/boundary/outer for each cell Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),& & gridctype(i,j)) End Do End do !$OMP end parallel do End subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! For each spline i determines if its support is inside, partially inside or outside of the simulation domain subroutine classifyspline(spl2,gridctype,bsptype) type(spline2d):: spl2 Integer:: gridctype(:,:) Integer:: bsptype(:) Integer:: nrank(2), nrz(2), ndegree(2) Integer:: i, j, mu, imin, imax, jmin, jmax Integer, Allocatable:: splinespan(:,:) call get_dim(spl2, nrank, nrz, ndegree) ! by default, all splines have part of their support outside the domain D bsptype=0 Allocate(splinespan(ndegree(1)+1,ndegree(2)+1)) Do mu=1,nrank(1)*nrank(2) ! scan the spline space ! obtain the axial spline index i=mod(mu-1,nrank(1))+1 ! obtain the radial spline index j=(mu-1)/nrank(1)+1 ! by default all cells are outside for correct behavior in boundaries splinespan=-1 ! find the axial span of this spline in cell indices imin=max(1,i-ndegree(1)) imax=min(nrz(1),i) ! find the radial span of this spline in cell indices jmin=max(1,j-ndegree(2)) jmax=min(nrz(2),j) ! obtain the cell type on which the spline is defined splinespan(1:(imax-imin+1),1:(jmax-jmin+1))=gridctype(imin:imax,jmin:jmax) if ( ANY( splinespan==1 ) ) then ! if at least one cell is fully in the domain the spline is an inside spline bsptype(mu)=1 else if (.not. ANY( splinespan==0 )) then ! if all the cells are outside, this is an outside spline bsptype(mu)=-1 end if end do End subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> general placeholder function used in fields_mod to compute the weight for weighted_bsplines !> This functions call the correct subroutine pointed by dirichlet_weight ! !--------------------------------------------------------------------------- SUBROUTINE geom_weight0(z,r,w) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w Real(kind=db):: rtmp(1:1), ztmp(1:1), wtmp(1:1,1:1) ztmp=z rtmp=r call Dirichlet_weight(ztmp,rtmp,wtmp) w=wtmp(1,1) End SUBROUTINE geom_weight0 SUBROUTINE geom_weight1(z,r,w) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w(0:) Real(kind=db):: rtmp(1:1), ztmp(1:1) Real(kind=db):: wtmp(0:size(w,1)-1,1:1) ztmp=z rtmp=r call Dirichlet_weight(ztmp,rtmp,wtmp) w=wtmp(:,1) End SUBROUTINE geom_weight1 SUBROUTINE geom_weight2(z,r,w) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) call Dirichlet_weight(z,r,w) ! End SUBROUTINE geom_weight2 !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> general placeholder function used in fields_mod to compute the domain weight for weighted_bsplines !> This functions call the correct subroutine pointed by domain_weight !> if the weight is negative return the id of the closest wall. !--------------------------------------------------------------------------- SUBROUTINE dom_weight0(z,r,w,idwall) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w Real(kind=db):: rtmp(1:1), ztmp(1:1), wtmp(1:1,1:1) INTEGER:: idwalltmp(1:1) INTEGER, optional, INTENT(OUT):: idwall ztmp=z rtmp=r call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp) w=wtmp(1,1) if (present(idwall)) idwall=idwalltmp(1) End SUBROUTINE dom_weight0 SUBROUTINE dom_weight1(z,r,w,idwall) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w(0:) Real(kind=db):: rtmp(1:1), ztmp(1:1) Real(kind=db):: wtmp(0:size(w,1)-1,1:1) INTEGER:: idwalltmp(1:1) INTEGER, optional, INTENT(OUT):: idwall ztmp=z rtmp=r call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp) w=wtmp(:,1) if (present(idwall)) idwall=idwalltmp(1) End SUBROUTINE dom_weight1 SUBROUTINE dom_weight2(z,r,w,idwall) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) INTEGER, optional, INTENT(OUT):: idwall(:) call domain_weight(z,r,w,idwall=idwall) ! End SUBROUTINE dom_weight2 SUBROUTINE dom_weight3(z,r,w,idwall) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:) Real(kind=db):: wtmp(1:1,size(w,1)) INTEGER, optional, INTENT(OUT):: idwall(:) call domain_weight(z,r,wtmp,idwall=idwall) !Write(*,*) 'wtmp, size', wtmp, size(wtmp) w=wtmp(1,:) ! End SUBROUTINE dom_weight3 SUBROUTINE geom_weightstd(z,r,w,wupper) ! return the geometric weight for a coaxial configuration ! or central conductor + ellipse if walltype ==1 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), OPTIONAL:: wupper(0:,:) Real(kind=db):: walltmp(0:size(w,1)-1,size(w,2)), elliptmp(0:size(w,1)-1,size(w,2)) Real(kind=db):: squareroot(size(w,2)), denom(size(w,2)) call cyllweight(z,r,walltmp, r_a, above1) SELECT CASE (walltype) CASE (0) call cyllweight(z,r,elliptmp, r_b, above2) CASE DEFAULT call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior) END SELECT if (present(wupper))then w=walltmp wupper=elliptmp return end if if(interior.eq.0 .and. above2.eq.0) then w=walltmp(:,:) else if (above1 .eq. 0) then w=elliptmp(:,:) else denom=walltmp(0,:)**2+elliptmp(0,:)**2 squareroot=sqrt(denom) w(0,:)=walltmp(0,:)+elliptmp(0,:)-squareroot ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(1,:)=walltmp(1,:)+elliptmp(1,:)-(elliptmp(1,:)*elliptmp(0,:)+walltmp(1,:)*walltmp(0,:))/& & squareroot ! z derivative of w w(2,:)=walltmp(2,:)+elliptmp(2,:)-(elliptmp(2,:)*elliptmp(0,:)+walltmp(2,:)*walltmp(0,:))/& & squareroot ! r derivative of w End If End if End SUBROUTINE geom_weightstd SUBROUTINE classification(z, r, ctype) ! classify if cell is fully inside, outside or on the boundary of the domain ! by calculating the weight on each corner and the cell. ! It is assumed that the cells are sufficiently small such that there is no sharp edge entering the cell ! where an inner portion of only one cell edge is outside of the domain real(kind=db), INTENT(IN):: r(2), z(2) INTEGER, INTENT(OUT):: ctype Real(kind=db)::weights(1:1,5) Real(kind=db):: zeval(5),reval(5) zeval=(/ z(1),z(2),z(1),z(2), (z(2)+z(1))/2 /) reval=(/ r(1),r(1),r(2),r(2), (r(2)+r(1))/2 /) CAll dom_weight(zeval,reval,weights) ctype=int(sign(1.0_db,weights(1,5))) If(weights(1,1)*weights(1,2) .le. 0 ) then ctype=0 return End If If(weights(1,1)*weights(1,3) .le. 0 ) then ctype=0 return End If If(weights(1,2)*weights(1,4) .le. 0 ) then ctype=0 return End If If(weights(1,3)*weights(1,4) .le. 0 ) then ctype=0 return End If end subroutine ! ################################################################## Subroutine calc_gauss(spl2, ngauss, i, j, xgauss, wgauss, gausssize, celltype) ! calculates the gauss integration points for the FEM method for any cell type ! takes care of boundary cells as well and limit the integration boundaries accordingly type(spline2d), INTENT(IN):: spl2 Integer, Intent(out):: gausssize INTEGER, Intent(out), Optional :: celltype Real(kind=db), Allocatable::rgrid(:),zgrid(:) - Real(kind=db), ALLOCATABLE::xgauss(:,:), wgauss(:) + Real(kind=db), ALLOCATABLE, intent(out)::xgauss(:,:), wgauss(:) Integer:: i,j, ngauss(2) Real(kind=db),Allocatable:: zpoints(:) Real(kind=db):: zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2)) Integer:: k, l, directiondown, directionup, nbzpoints, direction, ctype Logical:: hasmaxpoint Real(kind=db):: xptup, xptdown,wlu,wld, rmin, rmax type(spline1d):: splz, splr splz=spl2%sp1 splr=spl2%sp2 Allocate(zgrid(1:splz%nints+1)) Allocate(rgrid(1:splr%nints+1)) zgrid=splz%knots(0:splz%nints) rgrid=splr%knots(0:splr%nints) hasmaxpoint=.false. If(allocated(xgauss)) deallocate(xgauss) if(allocated(wgauss)) deallocate(wgauss) !Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),ctype) ctype=gridcelltype(i,j) If (ctype .ge. 1) then ! we have a normal internal cell Allocate(xgauss(ngauss(1)*ngauss(2),2)) Allocate(wgauss(ngauss(1)*ngauss(2))) gausssize=ngauss(1)*ngauss(2) !Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO Else If(ctype.eq.0) then ! we have a boundary cell directiondown=1 directionup=1 ! We check if the boundary goes through the cell upper and lower limit Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown) Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup) call dom_weight(zgrid(i),rgrid(j),wld) call dom_weight(zgrid(i),rgrid(j+1),wlu) select case ( directionup+directiondown) Case (0) ! The intersections are only on the left and right cell boundaries ! or the upper and lower limits are full boundaries nbzpoints=2 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),zgrid(i+1)/) Case(1) if(directiondown.eq.1)then if( wlu .gt. 0) then ! the lower left corner is inside nbzpoints=3 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i), xptdown,zgrid(i+1)/) else nbzpoints=2 Allocate(zpoints(nbzpoints)) if(wld.gt.0)then ! the upper left corner is inside zpoints=(/zgrid(i),xptdown/) else zpoints=(/xptdown,zgrid(i+1)/) end if end if else if(wld .gt. 0) then ! the lower left corner is inside nbzpoints=3 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i), xptup,zgrid(i+1)/) else nbzpoints=2 Allocate(zpoints(nbzpoints)) if(wlu.gt.0)then ! the upper left corner is inside zpoints=(/zgrid(i),xptup/) else zpoints=(/xptup,zgrid(i+1)/) end if end if end if Case(2) nbzpoints=4 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /) !If(wld.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /) !else ! nbzpoints=2 ! Allocate(zpoints(nbzpoints)) ! If(wld.ge.0) zpoints=(/zgrid(i),min(xptdown,xptup) /) ! If(wld.lt.0) zpoints=(/max(xptdown,xptup), zgrid(i+1) /) !end if End select Allocate(xgauss(ngauss(1)*ngauss(2)*(nbzpoints-1),2)) Allocate(wgauss(ngauss(1)*ngauss(2)*(nbzpoints-1))) gausssize=ngauss(1)*ngauss(2)*(nbzpoints-1) ! Compute gauss points Do l=1,nbzpoints-1 !CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) Call gauleg(zpoints(l),zpoints(l+1),zg,wzg,ngauss(1)) ! We test if the lower or upper side is in the domain call dom_weight(zg(1),rgrid(j),wld) rmin=rgrid(j) if (wld .le. 0) rmin = rgrid(j+1) Do k=1,ngauss(1) direction=2 Call Find_crosspointdico((/rgrid(j),rgrid(j+1)/),zg(k),rmax,direction) ! We compute the radial limits at each z position Call gauleg(min(rmin,rmax),max(rmin,rmax),rg,wrg,ngauss(2)) ! We obtain the gauss w and pos for these boundaries if(direction .eq. 0.and. wld .lt. 0) then wrg=0 end if xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),1) = zg(k) xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),2) = rg wgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1)) = wrg*wzg(k) End do End Do Else gausssize=1 Allocate(xgauss(1:1,2)) Allocate(wgauss(1:1)) !Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration xgauss(1,1)=(zg(1)+zg(ngauss(1)))*0.5 xgauss(1,2)=(rg(1)+rg(ngauss(2)))*0.5 wgauss(1)=0 End If If(PRESENT(celltype)) celltype=ctype End Subroutine calc_gauss Subroutine Find_crosspoint(x,y,xpt, direction) ! calculates the boundary limit between x(1) and x(2) using Newton's method Real(kind=db):: x(2), y Real(kind=db):: xptm, xpt, temp Real(kind=db):: w, wold Integer, Intent(INOUT):: direction Integer:: i ! calculates the position of the boundary ( where the weight changes sign ) ! between x(1) and x(2) at 2nd coordinate y xptm=x(1) xpt=x(2) i=0 ! direction=1 finds cross-point along z if(direction .eq.1) Call dom_weight(xptm,y,wold) ! direction=2 finds cross-point along r if(direction .eq.2) Call dom_weight(y,xptm,wold) if(direction .eq.1) Call dom_weight(xpt,y,w) if(direction .eq.2) Call dom_weight(y,xpt,w) ! if the weight doesn't change sign there is no cross point if(w*wold.gt.0) then direction=0 return End If ! Find the cross-point Do while( i .le.100 .and. abs(w).gt.1e-9) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold xptm=xpt-w*(xpt-xptm)/(w-wold) wold=w temp=xptm xptm=xpt xpt=temp i=i+1 if(direction .eq.1) Call dom_weight(xpt,y,w) if(direction .eq.2) Call dom_weight(y,xpt,w) End do if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0 End Subroutine Subroutine Find_crosspointdico(x,y,xpt, direction) ! calculates the boundary limit between x(1) and x(2) using dichotomy method Real(kind=db):: x(2), y Real(kind=db):: xpt Real(kind=db):: w1, w2, w3 Real(kind=db):: x1, x2, x3 Integer, Intent(INOUT):: direction Integer:: i ! calculates the position of the boundary ( where the weight changes sign ) ! between x(1) and x(2) at 2nd coordinate y x1=x(1) x2=x(2) i=0 select case(direction) case(1) ! direction=1 finds cross-point along z Call dom_weight(x1,y,w1) Call dom_weight(x2,y,w2) case(2) ! direction=2 finds cross-point along r Call dom_weight(y,x1,w1) Call dom_weight(y,x2,w2) end select ! if the weight doesn't change sign there is no cross point if(w1*w2.gt.0) then direction=0 xpt=x2 return End If ! Find the cross-point Do while( i .le. 500 .and. abs(x2-x1).gt.1e-13) x3=0.5*(x1+x2) i=i+1 select case(direction) case(1) ! direction=1 finds cross-point along z Call dom_weight(x3,y,w3) case(2) ! direction=2 finds cross-point along r Call dom_weight(y,x3,w3) end select if(w1*w3.gt.0)then! we are in a region were there is no change of sign x1=x3 w1=w3 else x2=x3 w2=w3 end if !if(abs(w3).lt.1e-14.and. w3 .ge.0)Exit End do xpt=x3 if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0 End Subroutine SUBROUTINE geom_rvaschtot(z,r,w,idwall) ! returns the total weight which is the same as the Dirichlet weight for a boundary ! defined with Rvaschev functions Use splinebound, ONLY: spline_w Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), allocatable:: w2(:,:) Real(kind=db), allocatable:: w3(:,:) INTEGER, optional, INTENT(OUT):: idwall(:) INTEGER:: sw2 sw2=size(w,2) if(present(idwall)) then idwall=0 allocate(w2(1:size(w,1),1:size(w,2))) allocate(w3(1:size(w,1),1:size(w,2))) call Dirichlet_weight(z,r,w2,w3) ! gives total weight where(w2(1,:).le.0) idwall=1 where(w3(1,:).le.0) idwall=2 call Combine(w2, w3, w, -1) else call Dirichlet_weight(z,r,w) end if End SUBROUTINE geom_rvaschtot Subroutine gstd(z,r,gtilde,w) ! g tilde function added to rhs of poisson solver for the standard coaxial configuration ! for the default weight function Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: belowtmp(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,1)) Real(kind=db):: denom(size(r,1)) if (above1.eq.0) then gtilde=0 gtilde(0,:)=Phiup RETURN end if ! Weight functions necessary for calculation of g ! coaxial insert call cyllweight(z,r,belowtmp, r_a, above1) SELECT CASE (walltype) CASE (0) ! top cylinder call cyllweight(z,r,abovetmp, r_b, above2) CASE DEFAULT ! Ellipse as in gt170 call ellipseweight(z,r,abovetmp,r_0,z_0,invr_z,invr_r,Interior) END SELECT ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(0,:)=(Phidown*abovetmp(0,:) + Phiup*belowtmp(0,:) ) / & & (abovetmp(0,:)+belowtmp(0,:)) If(size(gtilde,1) .gt. 2) then ! first derivative denom=(abovetmp(0,:)+belowtmp(0,:))**2 gtilde(1,:)=(Phiup-Phidown)*(belowtmp(1,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(1,:)) / & & denom ! Axial derivative gtilde(2,:)=(Phiup-Phidown)*(belowtmp(2,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(2,:)) / & & denom ! Radial derivative End If End subroutine SUBROUTINE gUpDown(z,r,gtilde,w) ! g tilde function added to rhs of poisson solver by combining two boundaries set at different potentials Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: belowtmp(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,1)) Real(kind=db):: denom(size(r,1)) ! Weight functions necessary for calculation of g call Dirichlet_weight(z,r,belowtmp,abovetmp) ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(0,:)=(Phidown*abovetmp(0,:) + Phiup*belowtmp(0,:) ) / & & (abovetmp(0,:)+belowtmp(0,:)) If(size(gtilde,2) .gt. 2) then ! first derivative denom=(abovetmp(0,:)+belowtmp(0,:))**2 gtilde(1,:)=(Phiup-Phidown)*(belowtmp(1,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(1,:)) / & & denom ! Axial derivative gtilde(2,:)=(Phiup-Phidown)*(belowtmp(2,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(2,:)) / & & denom ! Radial derivative End If END SUBROUTINE gUpDown Subroutine gtest(z,r,gtilde,w) ! calculates the Poisson gtilde term for testing the solver on a new geometry ! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: wtmp(0:size(gtilde,1)-1,size(gtilde,2)) ! !gtilde=0 !return if(present(w)) then wtmp=w else call Dirichlet_weight(z,r,wtmp) end if gtilde(0,:)=(sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))+2) If(size(gtilde,2) .gt. 1) then ! first derivative gtilde(1,:)=pi/(test_pars%Lz)*cos(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) gtilde(2,:)=pi/(test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*cos(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) End If gtilde(0,:)=gtilde(0,:)*(1-wtmp(0,:)) End subroutine Subroutine ftest(z,r,f) ! calculates the Poisson source term for testing the solver on a new geometry ! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT)::f(0:,:) !f(0,:)=1 !return ! f(0,:)=(pi/test_pars%Lz)**2*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))& & + (pi/test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))& & *( -1/r*(cos(pi*(r-test_pars%r0)/(test_pars%Lr))) + (pi/test_pars%Lr)*sin(pi*(r-test_pars%r0)/(test_pars%Lr))) End Subroutine ftest subroutine Reducematrix(full,reduced) ! Reduce the Finite element matrix in the LHS from full spline space to reduced web-spline space ! If A is the standard FEM matrix on the full grid space then the web-spline FEM matrix is etilde A etildet ! with etildet the transposed etilde matrix use mumps_bsplines type(mumps_mat):: full, reduced, tempmat1, tempmat2 Integer:: fullrank, reducedrank Integer:: i,j, k call to_mat(full) fullrank=full%rank reducedrank=nbreducedspline !WRITE(*,*) "# web-splines", nbreducedspline !WRITE(*,*) "# splines", fullrank call init(reducedrank,2,reduced) call init(fullrank,2,tempmat1) call init(reducedrank,2,tempmat2) tempmat1=mmx_mumps_mat_loc(full,etildet,(/fullrank,fullrank/),(/fullrank,reducedrank/),.false.) tempmat2=mmx_mumps_mat_loc(etildet,tempmat1,(/fullrank,reducedrank/),(/fullrank,reducedrank/),.true.) do i=1,reducedrank do k=tempmat2%irow(i),tempmat2%irow(i+1)-1 j=tempmat2%cols(k) call putele(reduced, i, j, tempmat2%val(k)) end do end do call to_mat(reduced) end subroutine !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION mmx_mumps_mat_loc(mata, matb, ranka, rankb, transpose) RESULT(matc) ! ! Return product matc=mata*matb ! All matrices are represented in csr sparse representation ! Use mumps_bsplines TYPE(mumps_mat) :: mata,matb,matc INTEGER:: ranka(2), rankb(2) ! nb of (rows,columns) for each matrix INTEGER :: n, info, m, request, sort LOGICAL:: transpose Character(1):: T ! m=ranka(1) n=ranka(2) T='N' if (transpose)then T='T' m=ranka(2) n=ranka(1) end if !#ifdef MKL ! if(associated(matc%val)) deallocate(matc%val) if(associated(matc%cols)) deallocate(matc%cols) if(associated(matc%irow)) deallocate(matc%irow) allocate(matc%irow(m+1)) allocate(matc%cols(1)) allocate(matc%val(1)) request=1 sort=7 CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow(1), & & matb%val, matb%cols, matb%irow(1), & & matc%val, matc%cols, matc%irow(1), & & 1, info) if(info .gt. 0) WRITE(*,'(a,i4)') " Error in mmx_mumps_mat_loc: ", info if(associated(matc%val)) deallocate(matc%val) if(associated(matc%cols)) deallocate(matc%cols) allocate(matc%val(matc%irow(m+1)-1)) allocate(matc%cols(matc%irow(m+1)-1)) request=2 CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow(1), & & matb%val, matb%cols, matb%irow(1), & & matc%val, matc%cols, matc%irow(1), & & 1, info) if(info .ne. 0) WRITE(*,'(a,i4)') " Error in mmx_mumps_mat_loc: ", info if (transpose)then matc%rank=ranka(2) else matc%rank=ranka(1) end if ! END FUNCTION mmx_mumps_mat_loc END MODULE geometry diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index c375154..3f217ca 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,938 +1,999 @@ MODULE splinebound USE constants USE bsplines USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints Use forSISLdata IMPLICIT NONE INTEGER, PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2 + logical:: nlexact type cellkind integer:: spldirkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: spltotkind=0 !< -1 outside (return -1) no dist to calculate; 0 boundary calculate dist with linked boundaries; 1 inside (return 1) no dist to calculate integer:: linkedboundaries(2)=0 !< stores the spline curve indices in the spline_domain of the spline boundaries that are the closest and at a distance lower than dist_extent (1) !< (1) is for dirichlet boundaries !< (2) is for domain boundaries integer:: leftknot(4)=0 !< knots pointer for s1424 in wtot then wdir real(kind=db):: lguess(2)=-1 !< Spline parameter left limit as start guess real(kind=db):: rguess(2)=-1 !< Spline parameter right limit as start guess end type cellkind TYPE spline_boundary ! all curves assume right handedness to set which side of the curve is inside or outside type(SISLCurve):: curve Real(kind=db):: Dirichlet_val !< Value for the dirichlet boundary condition created by this boundary Real(kind=db):: epsge=1.0e-5 !< geometric resolution used for calculating distances Real(kind=db):: epsce=1.0e-9 !< value of weight below which it is 0 INTEGER(kind(bd)):: type=bd_Dirichletconst !< type of boundary conditions END TYPE spline_boundary type spline_domain integer:: nbsplines = 0 !< number of spline boundaries in the domain type(spline_boundary), allocatable:: boundaries(:) !< List of boundaries in the domain Real(kind=db):: dist_extent=0.1 !< distance used for the merging with the plateau function for the weight type(cellkind), ALLOCATABLE:: cellk(:,:) !< Precomputed parameters at each cell for faster weight computation type(spline2d), pointer:: splrz => null() !< Pointer to the main spline grid used for the FEM solver Integer:: nb1 !< Number of grid points in the 1st dimension Integer:: nb2 !< Number of grid points in the 2nd dimension real(kind=db), ALLOCATABLE:: x1(:) !< Grid points in first direction for weight interpolation real(kind=db), ALLOCATABLE:: x2(:) !< Grid points in 2nd direction for weight interpolation real(kind=db), ALLOCATABLE:: dx1(:) !< inverse cell width in first direction for weight interpolation real(kind=db), ALLOCATABLE:: dx2(:) !< inverse cell width in 2nd direction for weight interpolation !type(SISLsurf):: Dirdomweight !< structure storing precalculated geometric weight for faster evaluation !type(SISLsurf):: totdomweight !< structure storing precalculated total weight for faster evaluation type(spline2d):: Dirdomweightspl !< structure storing precalculated geometric weight for faster evaluation type(spline2d):: totdomweightspl !< structure storing precalculated total weight for faster evaluation end type spline_domain CONTAINS !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Reads a spline domain from the namelist or from a h5 file. Needs to be called for the initialization of the module !> @param[in] Fileid Text file id of the input file containing namelists !> @param[out] spldom spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine read_splinebound(Fileid, spldom, splrz, rnorm, Phinorm) use mpi Integer:: Fileid type(spline_domain):: spldom type(spline2d):: splrz real(kind=db):: rnorm, phinorm Integer:: nbsplines, istat, mpirank, ierr real(kind=db):: dist_extent Character(len=128):: h5fname="", line real(kind=db) :: Dvals(30)=0 + integer:: i - namelist /spldomain/ nbsplines, dist_extent, h5fname, Dvals + namelist /spldomain/ nbsplines, dist_extent, h5fname, Dvals, nlexact CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) REWIND(fileid) READ(fileid,spldomain, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(fileid) read(fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, spldomain) Dvals=Dvals/phinorm dist_extent=dist_extent/rnorm if (.not. trim(h5fname)=='' ) then call setspline_domain(spldom, splrz, dist_extent, 0) call splinebound_readh5domain(h5fname,spldom, rnorm, phinorm) call classifycells(spldom) do i=1,spldom%nbsplines spldom%boundaries(i)%Dirichlet_val=Dvals(i) end do return else WRITE(*,*) "Error the filename h5fname is not defined. No boundary has been set!" call mpi_Abort(MPI_COMM_WORLD, -1, ierr) end if end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Saves the spline boundaries to the result file !> @param[in] File_handle futils h5 file id !> @param[in] curr_grp groupname under which the boundaries must be saved !> @param[in] spldom spline domain !--------------------------------------------------------------------------- Subroutine splinebound_diag(File_handle, curr_grp, spldom) use mpi Use futils Use basic, ONLY: rnorm, phinorm Integer:: File_handle type(spline_domain):: spldom Character(len=*):: curr_grp CHARACTER(len=128):: grpname Integer:: ierr, mpirank, i CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0) THEN Write(grpname,'(a,a)') trim(curr_grp),"/geometry_spl" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "dist_extent",spldom%dist_extent) Call attach(File_handle, trim(grpname), "nbsplines", spldom%nbsplines) do i=1,spldom%nbsplines Write(grpname,'(a,a,i2.2)') trim(curr_grp),"/geometry_spl/",i If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "Dirichlet_val", spldom%boundaries(i)%Dirichlet_val*phinorm) Call attach(File_handle, trim(grpname), "order", spldom%boundaries(i)%curve%ik) Call attach(File_handle, trim(grpname), "kind", spldom%boundaries(i)%curve%ikind) Call attach(File_handle, trim(grpname), "type", spldom%boundaries(i)%type) Call attach(File_handle, trim(grpname), "dim", spldom%boundaries(i)%curve%idim) CALL putarr(File_handle, TRIM(grpname)//"/pos", spldom%boundaries(i)%curve%ecoef*rnorm) CALL putarr(File_handle, TRIM(grpname)//"/knots", spldom%boundaries(i)%curve%et) end do END IF End subroutine splinebound_diag !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Read a spline boundary domain from an h5 file structure !> @param[out] spldom new spline domain !> @param[in] filename filename of the h5 file !> @param[in] rnorm distance normalization constant !> @param[in] phinorm electric potential normalization constant !--------------------------------------------------------------------------- subroutine splinebound_readh5domain(filename, spldom, rnorm, phinorm) use futils use forSISL implicit none Character(len=*),intent(in) :: filename type(spline_domain),intent(inout) :: spldom integer:: h5id, i real(kind=db):: rnorm, phinorm CHARACTER(len=128):: grpname integer:: periodic integer:: order, dim, bdtype INTEGER:: posrank, posdim(2), err real(kind=db):: Dval, epsge, epsce real(kind=db),allocatable:: points(:,:) call openf(filename, h5id,'r','d') call getatt(h5id, '/geometry_spl/','nbsplines', spldom%nbsplines) ! prepare memory if (allocated(spldom%boundaries)) then do i=1,size(spldom%boundaries,1) call free_bsplinecurve(spldom%boundaries(i)) end do DEALLOCATE(spldom%boundaries) end if allocate(spldom%boundaries(spldom%nbsplines)) ! Read each boundary curve individually do i=1,spldom%nbsplines Write(grpname,'(a,i2.2)') "/geometry_spl/",i If(.not. isgroup(h5id, trim(grpname))) THEN Write(*,*) "Error the geometry definition file is invalid" END IF periodic=0 Call getatt(h5id, trim(grpname), "Dirichlet_val", Dval) Call getatt(h5id, trim(grpname), "epsge", epsge) Call getatt(h5id, trim(grpname), "epsce", epsce) Call getatt(h5id, trim(grpname), "order", order) Call getatt(h5id, trim(grpname), "dim", dim) err=0 Call getatt(h5id, trim(grpname), "periodic", periodic,err) if(err .lt.0) periodic=0 CALL getdims(h5id, TRIM(grpname)//"/pos", posrank, posdim) allocate(points(posdim(1),posdim(2))) CALL getarr(h5id, TRIM(grpname)//"/pos", points) points=points/rnorm + Call setspline_boundary(spldom%boundaries(i),transpose(points), order-1, Dval/phinorm, epsge,epsce, periodic) bdtype=bd err=0 Call getatt(h5id, trim(grpname), "type", bdtype,err) if(err.ge.0) spldom%boundaries(i)%type=bdtype deallocate(points) end do call closef(h5id) end subroutine splinebound_readh5domain !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> initialize a spline domain and allocate the necessary memory !> @param[out] spldom new spline domain !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] dist_extent normalized characteristic fall lenght of the weight !> @param[in] nb_splines number of boundary splines to allocate !--------------------------------------------------------------------------- subroutine setspline_domain(spldom,splrz,dist_extent, nb_splines) type(spline_domain):: spldom type(spline2d), TARGET:: splrz real(kind=db):: dist_extent integer:: nb_splines, nb1, nb2 ! Store the grid parameters to speed-up calculations nb1=splrz%sp1%nints nb2=splrz%sp2%nints spldom%nb1=nb1 spldom%nb2=nb2 spldom%splrz=>splrz allocate(spldom%cellk(0:nb1-1,0:nb2-1)) allocate(spldom%x1(0:nb1)) allocate(spldom%x2(0:nb2)) allocate(spldom%dx1(0:nb1-1)) allocate(spldom%dx2(0:nb2-1)) spldom%x1(0:)=splrz%sp1%knots(0:nb1) spldom%x2(0:)=splrz%sp2%knots(0:nb2) spldom%dx1(0:)=1/(spldom%x1(1:nb1)-spldom%x1(0:nb1-1)) spldom%dx2(0:)=1/(spldom%x2(1:nb2)-spldom%x2(0:nb2-1)) !Prepare structures to host singular spline boundaries spldom%nbsplines=nb_splines if(spldom%nbsplines.gt. 0) allocate(spldom%boundaries(nb_splines)) spldom%dist_extent=dist_extent end subroutine setspline_domain !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief initialize a spline boundary and allocate the necessary memory !> @param[out] b_curve new spline boundary !> @param[in] cpoints control points at the node positions !> @param[in] degree degree of the spline polynomia defining the boundary curve !> @param[in] D_val Normalized value of the Dirichlet boundary condition for this curve !> @param[in] epsge geometric precision used by SISL !> @param[in] epsce arithmetic precision used by SISL !> @param[in] periodic set if the spline curve is periodic !--------------------------------------------------------------------------- subroutine setspline_boundary(b_curve, cpoints, degree, D_val, epsge, epsce, periodic) Use bsplines use forSISL,ONLY: newcurve, s1630 use mpi type(spline_boundary):: b_curve Real(kind=db):: cpoints(:,:) Real(REAL64),ALLOCATABLE:: points(:) Real(REAL64):: astpar integer:: degree integer, optional:: periodic - Integer:: order, ierr - Real(kind=db):: D_val + Integer:: order, ierr, i,j + Real(kind=db):: D_val ,dist Real(kind=db),OPTIONAL :: epsge, epsce Integer:: nbpoints, dim, jstat, bsptype integer:: period period=0 if(present(periodic))period=periodic nbpoints= size(cpoints,2) dim=size(cpoints,1) order=degree+1 if(nbpoints .lt. order) then WRITE(*,'(a,i3,a,i5)') "Error: the number of points", nbpoints, " is insuficient for the required order ", order CALL mpi_finalize(ierr) call EXIT(-1) end if allocate(points(dim*nbpoints)) - points=reshape(cpoints,(/dim*nbpoints/)) + j=1 + points(1:2)=cpoints(:,1) + do i=2,nbpoints + dist=sum((points(2*(j-1)+1:2*(j-1)+2)-cpoints(:,i))**2) + if(dist.lt.1e-12) cycle + points(2*j+1:2*j+2)=cpoints(:,i) + j=j+1 + end do + + + !points=reshape(cpoints,(/dim*nbpoints/)) bsptype=1 ! open boundaries b-spline if(period.gt.0) bsptype=-1 ! closed periodic curve - astpar=0.0 ! starting parameter for the knots vector + astpar=0.0_db ! starting parameter for the knots vector ! initialize a new curve using SISL - CALL s1630(points, nbpoints, astpar, bsptype, dim, order, b_curve%curve, jstat) + CALL s1630(points, j, astpar, bsptype, dim, order, b_curve%curve, jstat) if (jstat > 0 ) WRITE(*,*) "Warning ", jstat," in curve initialisation s1630 for splineweight" if (jstat < 0 ) WRITE(*,*) "Error ", jstat," in curve initialisation s1630 for splineweight" b_curve%Dirichlet_val=D_val if(present(epsge)) b_curve%epsge=epsge if(present(epsce)) b_curve%epsce=epsce end subroutine setspline_boundary !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the Dirichlet boundary weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_w(spldom,x1,x2,w) - use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(0:,:) Integer,allocatable::i(:),j(:) + Integer:: k,l + allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) + if (nlexact) then + do k=1,size(x1) + l=spldom%cellk(i(k),j(k))%linkedboundaries(1) + if(l.eq.0)then + w(:,k)=0 + w(0,k)=spldom%cellk(i(k),j(k))%spldirkind + else + call splineweight(spldom%boundaries(l),x1(k),x2(k),w(:,k),spldom%dist_extent) + end if + end do + return + end if + if (size(w,1).gt.1) then CALL speval(spldom%Dirdomweightspl, x1, x2, i, j, w(0,:), w(1,:), w(2,:)) else CALL speval(spldom%Dirdomweightspl, x1, x2, i, j, w(0,:)) end if End SUBROUTINE spline_w !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the total geometric weight from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] x1(:) array of axial positions where the weights are evaluated !> @param[in] x2(:) array of radial positions where the weights are evaluated !> @param[out] w(:,0:) matrix of weights with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_wtot(spldom,x1,x2,w,idwall) use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(0:,:) INTEGER, optional, INTENT(OUT):: idwall(:) - Integer:: k + Integer:: k,l Integer,allocatable::i(:),j(:) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) + if(present(idwall)) then Do k=1,size(x2,1) idwall(k)=spldom%cellk(i(k),j(k))%linkedboundaries(2) END DO end if + if (nlexact) then + do k=1,size(x1) + l=spldom%cellk(i(k),j(k))%linkedboundaries(2) + if(l.eq.0)then + w(:,k)=0 + w(0,k)=spldom%cellk(i(k),j(k))%spltotkind + else + call splineweight(spldom%boundaries(l),x1(k),x2(k),w(:,k),spldom%dist_extent) + end if + end do + return + end if + if (size(w,1).gt.1) then CALL speval(spldom%totdomweightspl, x1, x2, i, j, w(0,:), w(1,:), w(2,:)) else CALL speval(spldom%totdomweightspl, x1, x2, i, j, w(0,:)) end if End SUBROUTINE spline_wtot !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the interpolation in the domain of the Dirichlet boundary conditions from a given spline domain !> @param[in] spldom spline domain containing the information on the boundary conditions !> @param[in] z(:) array of axial positions where the weights are evaluated !> @param[in] r(:) array of radial positions where the weights are evaluated !> @param[out] g(:,0:) matrix of boundary interpolations g with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- SUBROUTINE spline_g(spldom,x1,x2,g,w) use forSISL,ONLY: s1424 use bsplines type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: g(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) REAL(real64),allocatable:: gtmp(:,:) - Integer:: k + Integer:: k,l Integer,allocatable::i(:),j(:) !type(cellkind):: cellk allocate(gtmp(0:size(g,1)-1,size(x2,1))) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if(present(w)) then gtmp=w else - CALL speval(spldom%Dirdomweightspl, x1, x2,i,j, gtmp(0,:), gtmp(1,:), gtmp(2,:)) + if (nlexact) then + do k=1,size(x1) + l=spldom%cellk(i(k),j(k))%linkedboundaries(2) + call splineweight(spldom%boundaries(l),x1(k),x2(k),gtmp(:,k),spldom%dist_extent) + end do + else + CALL speval(spldom%Dirdomweightspl, x1, x2,i,j, gtmp(0,:), gtmp(1,:), gtmp(2,:)) + end if end if Do k=1,size(x2,1) if(spldom%cellk(i(k),j(k))%spldirkind.eq.0)then if(gtmp(0,k) .ge. 0) then if(size(g,1) .gt. 1) then g(1:2,k)=-gtmp(1:2,k)*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val end if g(0,k)=(1-gtmp(0,k))*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val else g(0,k)=spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val if(size(g,1).gt. 1) then g(1:2,k)=0 end if end if else g(:,k)=0 end if end DO End SUBROUTINE spline_g !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Evaluates the geometric weight induced by the spline curve defined by b_curve at position (z,r) !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] z axial position where the weight is evaluated !> @param[in] r radial position where the weight is evaluated !> @param[out] weight(:) weight index defines the order of derivation by r or z !> @param[in] h distance from the spline at which the weight is 1 !> @param[out] distance unscaled distance between evaluation point and spline b_curve !> @param[inout] leftknot initial guess for the closest spline knot of the points (r,z) !--------------------------------------------------------------------------- subroutine splineweight(b_curve, z, r, weight, h, distance, guess, lguess, rguess) Use forSISL, ONLY: s1227,s1221, s1774 type(spline_boundary):: b_curve Real(kind=db)::r,z Real(kind=db):: weight(0:) Real(kind=db),OPTIONAL:: distance real(kind=db),OPTIONAL:: guess real(kind=db),OPTIONAL:: lguess real(kind=db),OPTIONAL:: rguess integer:: sstatus, der, left,siz real(kind=db):: h, d, tpos, proj, norm real(kind=real64):: curvepos(2*b_curve%curve%idim) real(kind=db):: leftpar, rightpar,guesspar weight=0 der=1 sstatus=-1 guesspar=-1.0_db if(present(lguess) .and. present(rguess)) then leftpar=lguess rightpar=rguess guesspar=(lguess+rguess)/2 call s1774(b_curve%curve,(/z,r/),b_curve%curve%idim,b_curve%epsge,leftpar,rightpar,guesspar,tpos,sstatus) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1774 for splineweight at ", z, r else call dist(b_curve,(/z,r/),d,tpos) end if ! position and derivative wrt r,z - call s1221(b_curve%curve,der,tpos,left,curvepos,sstatus) + call s1227(b_curve%curve,der,tpos,left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", z, r if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", z, r d=sqrt((curvepos(1)-z)**2+(curvepos(2)-r)**2) - weight(0)=1-max((h-d)/h,0.0_db)**3 + weight(0)=1-max((h-d)/h,0.0_db)**2 norm=sqrt(curvepos(3)**2+curvepos(4)**2) if(norm.gt.0) curvepos(3:4)=curvepos(3:4)/norm ! if the projection of the distance vector on the normal is negative, the weight is negative proj=(-(z-curvepos(1))*curvepos(4)+(r-curvepos(2))*curvepos(3)) - if (proj .lt. 0 .or. abs(abs(proj) -sqrt((z-curvepos(1))**2+(r-curvepos(2))**2)).gt.1e-10) weight(0)=-weight(0) + !if (proj .lt. 0 .or. abs(abs(proj) -sqrt((z-curvepos(1))**2+(r-curvepos(2))**2)).gt.1e-8) weight(0)=-weight(0) + if (proj .lt. 0)then + weight(0)=-weight(0) + end if !if (proj .lt. 0 ) weight(0)=-weight(0) siz=size(weight,1) if (size(weight,1).gt.1 .and. abs(weight(0)) .lt. 1) then - weight(1)=-3*curvepos(4)*(h-d)**2/h**3 - weight(2)=+3*curvepos(3)*(h-d)**2/h**3 + weight(1)=-2*curvepos(4)*abs((h-d))/h**2 + weight(2)=+2*curvepos(3)*abs((h-d))/h**2 end if if(present(distance)) distance=d if(present(guess)) guess=tpos end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the closest distance between the point and the selected spline b_curve !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] point(:) array containing the position from which to calculate the distance !> @param[out] distance distance from the point to the spline !> @param[in] pos parameter value of the closest point on the spline !--------------------------------------------------------------------------- subroutine dist(b_curve, point, distance, pos) Use forSISL, ONLY: s1957,s1953, s1221,s1227 type(spline_boundary):: b_curve Real(kind=db):: point(:) real(kind=db):: distance Real(kind=db),optional::pos REAL(real64):: posres, epsco, epsge,curvepos(2),d,distmin REAL(real64),allocatable::intpar(:) integer:: numintpt, numintcu,i,left,sstatus type(SISLIntCurve),ALLOCATABLE:: intcurve(:) epsco=1.0e-15 epsge=1.0e-15 !epsco=0 !epsge=b_curve%epsge numintpt=0 sstatus=0 distmin=HUGE(d) + !call s1957(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,posres,distance,sstatus) +! + !if(sstatus.eq.0) then + ! if (present(pos)) pos=posres + ! return + !end if +! + !!if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1953 for splineweight at ", point(1), point(2) + !if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1953 for splineweight at ",point(1), point(2) +! call s1953(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,numintpt,intpar,numintcu,intcurve,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1953 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1953 for splineweight at ",point(1), point(2) if(numintpt .gt. 1) then Do i=1,numintpt call s1227(b_curve%curve,0,intpar(i),left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1221 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1221 for splineweight at ",point(1), point(2) d=(curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2 if(d .lt. distmin) then distmin=d posres=intpar(i) end if end do else if(numintpt .gt. 0) then posres=intpar(1) end if distance=distmin if(numintcu.ge.1) then posres=0.5*(intcurve(1)%epar1(1)+intcurve(1)%epar1(2)) end if - call s1221(b_curve%curve,0,posres,left,curvepos,sstatus) + call s1227(b_curve%curve,0,posres,left,curvepos,sstatus) if (sstatus > 0 ) WRITE(*,*) "Warning ",sstatus," in distance calculation s1227 for splineweight at ", point(1), point(2) if (sstatus < 0 ) WRITE(*,*) "Error ",sstatus," in distance calculation s1227 for splineweight at ", point(1), point(2) distance=sqrt((curvepos(1)-point(1))**2+(curvepos(2)-point(2))**2) if (present(pos)) pos=posres END subroutine SUBROUTINE classify(x1, x2, cellk, spldom, wpredir, wpretot) real(kind=db), INTENT(IN):: x2(2), x1(2) type(cellkind), intent(INOUT):: cellk type(spline_domain)::spldom Real(kind=db):: zeval(4),reval(4), wpretot, wpredir real(kind=db), allocatable:: guess(:,:), w(:,:,:) Real(kind=db):: dmin, insidedir, insidetot, distance integer:: i,k allocate(guess(spldom%nbsplines,4)) allocate(w(0:2,spldom%nbsplines,4)) w=0 cellk%spldirkind=0 guess=-1.0_db dmin=HUGE(spldom%dist_extent) cellk%linkedboundaries=0 zeval=(/ x1(1),x1(2),x1(1),x1(2) /) reval=(/ x2(1),x2(1),x2(2),x2(2) /) insidedir=1 insidetot=1 do i=1,spldom%nbsplines do k=1,4 ! calculate the weight for each spline boundaries at each cell corner call splineweight(spldom%boundaries(i),zeval(k),reval(k),w(:,i,k),spldom%dist_extent,distance,guess(i,k)) ! We find the closest boundary to this point if(distance .lt. dmin) then ! If we are close enough we check if we are below dist_extent and need to calculate the distance each time if(distance .lt. spldom%dist_extent) then if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then cellk%linkedboundaries(1)=i cellk%spldirkind=0 end if cellk%linkedboundaries(2)=i cellk%spltotkind=0 end if dmin=distance ! Otherwise we define the interior by the closest spline if(spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) then insidedir=w(0,i,k) end if insidetot=w(0,i,k) end if ! The neumann boundaries take precedence over the dirichlet boundaries ! this is important when they define what is outside of the simulation domain. if(spldom%boundaries(i)%type.eq. bd_Neumann.and. w(0,i,k).lt.0)then insidetot=w(0,i,k) if(distance.lt.spldom%dist_extent)then cellk%linkedboundaries(2) =i else cellk%linkedboundaries(2) =0 end if end if end do end do if(cellk%linkedboundaries(1) .gt. 0) then i=cellk%linkedboundaries(1) cellk%lguess(1)=minval(guess(i,:),1,guess(i,:).ge.0) cellk%rguess(1)=maxval(guess(i,:),1) wpredir=w(0,i,1) else cellk%spldirkind=sign(1,int(insidedir)) wpredir=insidedir end if if(cellk%linkedboundaries(2) .gt. 0) then i=cellk%linkedboundaries(2) wpretot=w(0,i,1) else cellk%spltotkind=sign(1,int(insidetot)) wpretot=insidetot end if end subroutine subroutine classifycells(spldom) use forSISL, ONLY: s1537, s1424 use bsplines type(spline_domain):: spldom integer:: i,j, dims(2), nbeval1, nbeval2,k,l real(kind=db)::val type(cellkind):: cellk real(kind=db), allocatable:: wpretot(:,:,:), wpredir(:,:,:), c(:,:), x1(:), x2(:) allocate(wpretot(1:1,0:spldom%nb1,0:spldom%nb2)) allocate(wpredir(1:1,0:spldom%nb1,0:spldom%nb2)) nbeval1=spldom%nb1+3 nbeval2=spldom%nb2+3 ! We set the interpolation points such that the spline interpolation of the weight uses the same knots as the spline interpolation of the electric potential allocate(x1(0:nbeval1-1),x2(0:nbeval2-1)) x1(0)=spldom%x1(0) x1(1)=(spldom%x1(0)+spldom%x1(1))/2.0_db j=0 do i=2,spldom%nb1 j=j+1 x1(i)=spldom%x1(j) !x1(i)=2*spldom%x1(j)-x1(i-1) end do x1(nbeval1-2)=(spldom%x1(spldom%nb1-1)+3*spldom%x1(spldom%nb1))/2.0_db x1(nbeval1-1)=spldom%x1(spldom%nb1) !write(*,*)"x1", x1 ! We do the same for x2 x2(0)=spldom%x2(0) x2(1)=(spldom%x2(0)+spldom%x2(1))/2.0_db j=0 do i=2,spldom%nb2 j=j+1 x2(i)=spldom%x2(j) !x2(i)=2*spldom%x2(j)-x2(i-1) end do x2(nbeval2-2)=(spldom%x2(spldom%nb2-1)+spldom%x2(spldom%nb2))/2.0_db x2(nbeval2-1)=spldom%x2(spldom%nb2) !write(*,*)"x2", x2 wpretot=0 wpredir=0 !$OMP PARALLEL DO private(i,j) do i=0,spldom%nb1-1 !DIR$ UNROLL do j=0,spldom%nb2-1 call classify(spldom%x1(i:i+1),spldom%x2(j:j+1),spldom%cellk(i,j),spldom, wpredir(1,i,j),wpretot(1,i,j)) end do end do !$OMP END PARALLEL DO deallocate(wpretot) deallocate(wpredir) allocate(wpretot(1:1,0:nbeval1-1,0:nbeval2-1)) allocate(wpredir(1:1,0:nbeval1-1,0:nbeval2-1)) !$OMP PARALLEL DO private(i,j,cellk,k,l) do i=0,nbeval1-1 call locintv(spldom%splrz%sp1,x1(i),k) do j=0,nbeval2-1 call locintv(spldom%splrz%sp2,x2(j),l) cellk=spldom%cellk(k,l) If(abs(cellk%spldirkind) .eq. 1) Then wpredir(1,i,j)=cellk%spldirkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(i),x2(j), wpredir(:,i,j),spldom%dist_extent) end IF If(abs(cellk%spltotkind) .eq. 1) Then wpretot(1,i,j)=cellk%spltotkind else call splineweight(spldom%boundaries(cellk%linkedboundaries(2)), x1(i),x2(j), wpretot(:,i,j),spldom%dist_extent) end IF end do end do !$OMP END PARALLEL DO ! Set the approximated spline weight for the Dirichlet boundary conditions CALL set_splcoef((/3,3/),x1,x2,spldom%Dirdomweightspl) call get_dim(spldom%Dirdomweightspl,dims) !Write(*,*) "size x1, x2 knots", size(x1),size(x2),dims, size(wpredir) allocate(c(dims(1),dims(2))) call get_splcoef(spldom%Dirdomweightspl, wpredir(1,:,:), c) CALL gridval(spldom%Dirdomweightspl,spldom%x1(1),spldom%x2(1), val ,(/0,0/),c) !write(*,*)"x2", x2 !write(*,*)"konot1", spldom%x1 !write(*,*)"konots1 interp", spldom%Dirdomweightspl%sp1% ! Set the approximated spline weight for the Neumann boundary conditions CALL set_splcoef((/3,3/),x1,x2,spldom%totdomweightspl) call get_splcoef(spldom%totdomweightspl, wpretot(1,:,:), c) CALL gridval(spldom%totdomweightspl,spldom%x1(1),spldom%x2(1), val ,(/0,0/),c) deallocate(c) end subroutine subroutine getindex(x1,x2,spldom, i, j) use distrib, ONLY: closest type(spline_domain):: spldom real(kind=db):: x1(:), x2(:) integer:: i(:),j(:) call locintv(spldom%splrz%sp1,x1, i) call locintv(spldom%splrz%sp2,x2, j) end subroutine SUBROUTINE speval(sp, xp, yp, leftx, lefty, f00, f10, f01) ! ! Compute the function f00 and its derivatives ! f10 = d/dx f ! f01 = d/dy f ! assuming that its PPFORM/BCOEFSC was already computed! ! TYPE(spline2d), INTENT(inout) :: sp DOUBLE PRECISION, DIMENSION(:), INTENT(in) :: xp, yp INTEGER, DIMENSION(:), INTENT(in) :: leftx, lefty DOUBLE PRECISION, DIMENSION(:), INTENT(out) :: f00 DOUBLE PRECISION, DIMENSION(:), INTENT(out), OPTIONAL :: f10, f01 ! INTEGER :: np DOUBLE PRECISION :: x(SIZE(xp)), y(SIZE(yp)) INTEGER :: i, nidbas(2) DOUBLE PRECISION :: temp0(SIZE(xp),sp%sp2%order), temp1(SIZE(xp),sp%sp2%order) LOGICAL :: nlppform ! ! Apply periodicity if required ! np = SIZE(xp) nidbas(1) = sp%sp1%order-1 nidbas(2) = sp%sp2%order-1 nlppform = sp%sp1%nlppform .OR. sp%sp2%nlppform ! ! Locate the interval containing x, y ! x(:) = xp(:) - sp%sp1%knots(leftx(:)) y(:) = yp(:) - sp%sp2%knots(lefty(:)) ! ! Compute function/derivatives ! ! Using PPFORM !---------- DO i=1,np CALL my_ppval1(nidbas(1), x(i), sp%ppform(:,leftx(i)+1,:,lefty(i)+1), & & temp0(i,:), temp1(i,:)) END DO ! CALL my_ppval0(nidbas(2), y, temp0, 0, f00) if(present(f01))then CALL my_ppval0(nidbas(2), y, temp0, 1, f01) end if if(present(f10))then CALL my_ppval0(nidbas(2), y, temp1, 0, f10) end if !----------- CONTAINS !+++ SUBROUTINE my_ppval0(p, x, ppform, jder, f) ! ! Compute function and derivatives from the PP representation ! for many points x(:) INTEGER, INTENT(in) :: p DOUBLE PRECISION, INTENT(in) :: x(:) DOUBLE PRECISION, INTENT(in) :: ppform(:,:) INTEGER, INTENT(in) :: jder DOUBLE PRECISION, INTENT(out) :: f(:) DOUBLE PRECISION :: fact INTEGER :: j SELECT CASE (jder) CASE(0) ! function value SELECT CASE(p) CASE(1) f(:) = ppform(:,1) + x(:)*ppform(:,2) CASE(2) f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*ppform(:,3)) !!$ CASE(3) !!$ f(:) = ppform(:,1) + x(:)*(ppform(:,2)+x(:)*(ppform(:,3)+x(:)*ppform(:,4))) CASE(3:) f(:) = ppform(:,p+1) DO j=p,1,-1 f(:) = f(:)*x(:) + ppform(:,j) END DO END SELECT CASE(1) ! 1st derivative SELECT CASE(p) CASE(1) f(:) = ppform(:,2) CASE(2) f(:) = ppform(:,2) + x(:)*2.d0*ppform(:,3) !!$ CASE(3) !!$ f(:) = ppform(:,2) + x(:)*(2.d0*ppform(:,3)+x(:)*3.0d0*ppform(:,4)) CASE(3:) f(:) = p*ppform(:,p+1) DO j=p-1,1,-1 f(:) = f(:)*x(:) + j*ppform(:,j+1) END DO END SELECT CASE default ! 2nd and higher derivatives f(:) = ppform(:,p+1) fact = p-jder DO j=p,jder+1,-1 f(:) = f(:)/fact*j*x(:) + ppform(:,j) fact = fact-1.0d0 END DO DO j=2,jder f(:) = f(:)*j END DO END SELECT END SUBROUTINE my_ppval0 !+++ SUBROUTINE my_ppval1(p, x, ppform, f0, f1) ! ! Compute function and first derivative from the PP representation INTEGER, INTENT(in) :: p DOUBLE PRECISION, INTENT(in) :: x DOUBLE PRECISION, INTENT(in) :: ppform(:,:) DOUBLE PRECISION, INTENT(out) :: f0(:) DOUBLE PRECISION, INTENT(out) :: f1(:) DOUBLE PRECISION :: fact INTEGER :: j SELECT CASE(p) CASE(1) f0(:) = ppform(1,:) + x*ppform(2,:) f1(:) = ppform(2,:) CASE(2) f0(:) = ppform(1,:) + x*(ppform(2,:)+x*ppform(3,:)) f1(:) = ppform(2,:) + x*2.d0*ppform(3,:) CASE(3) f0(:) = ppform(1,:) + x*(ppform(2,:)+x*(ppform(3,:)+x*ppform(4,:))) f1(:) = ppform(2,:) + x*(2.d0*ppform(3,:)+x*3.0d0*ppform(4,:)) CASE(4:) f0 = ppform(p+1,:) f1 = f0 DO j=p,2,-1 f0(:) = ppform(j,:) + x*f0(:) f1(:) = f0(:) + x*f1(:) END DO f0(:) = ppform(1,:) + x*f0(:) END SELECT END SUBROUTINE my_ppval1 !+++ END SUBROUTINE speval subroutine free_bsplinecurve(b_curve) type(spline_boundary):: b_curve call freeCurve(b_curve%curve) !call freeIntCurve(b_curve%intcurve) end subroutine END MODULE splinebound