diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index 64d2bed..d568c93 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,1031 +1,1042 @@ MODULE splinebound USE constants USE bsplines USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints Use forSISLdata IMPLICIT NONE INTEGER(SELECTED_INT_KIND(4)), PARAMETER :: bd=-1, bd_Dirichletconst=0, bd_Dirichletvar=1, bd_Neumann=2 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, Potinn, Potout Character(len=128):: h5fname="", line real(kind=db) :: Dvals(30)=0 integer:: nbpoints integer:: degree=4, i namelist /spldomain/ nbsplines, dist_extent, h5fname, Dvals 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), "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 Real(kind=db),OPTIONAL :: epsge, epsce Integer:: nbpoints,i, 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/)) bsptype=1 ! open boundaries b-spline if(period.gt.0) bsptype=-1 ! closed periodic curve astpar=0.0 ! starting parameter for the knots vector ! initialize a new curve using SISL CALL s1630(points, nbpoints, 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:: k,sstat,l !type(cellkind):: cellk REAL(real64):: wtmp(4), point(2) Integer,allocatable::i(:),j(:) real(kind=db),allocatable:: wtmp2(:,:) allocate(i(size(x2,1)),j(size(x2,1))) allocate(wtmp2(size(x1,1),3)) call getindex(x1, x2, spldom, i, j) !Do k=1,size(x2,1) ! ! cellk=spldom%cellk(i(k),j(k)) ! ! SELECT CASE (cellk%spldirkind) ! CASE(-1) ! w(k,:)=0 ! w(k,0)=-1 ! CYCLE ! CASE(1) ! w(k,:)=0 ! w(k,0)=1 ! CYCLE ! CASE DEFAULT ! !call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1)) ! ! !point=(/x1(k),x2(k)/) ! !call s1424(spldom%Dirdomweight,1,1,point,cellk%leftknot(3),cellk%leftknot(4),wtmp,sstat) ! !if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight" ! !if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight" ! !w(k,:)=wtmp(1:size(w,2)) ! ! ! END SELECT !end DO if (size(w,2).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 !Do k=1,size(x2,1) ! SELECT CASE (spldom%cellk(i(k),j(k))%spldirkind) ! CASE(-1) ! w(k,:)=0 ! w(k,0)=-1 ! CYCLE ! CASE(1) ! w(k,:)=0 ! w(k,0)=1 ! CYCLE ! CASE DEFAULT ! end select !end do 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,sstat,l REAL(real64):: wtmp(4), point(2) Integer,allocatable::i(:),j(:),celltype(:) real(kind=db),allocatable:: wtmp2(:,:) 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 !Do k=1,size(x2,1) ! ! cellk=spldom%cellk(i(k),j(k)) ! if(present(idwall)) idwall(k)=cellk%linkedboundaries(2) ! SELECT CASE (cellk%spltotkind) ! CASE(-1) ! w(k,:)=0 ! w(k,0)=-1 ! CYCLE ! CASE(1) ! w(k,:)=0 ! w(k,0)=1 ! CYCLE ! CASE DEFAULT ! !call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,rguess=cellk%rguess(1),lguess=cellk%lguess(1)) ! ! point=(/x1(k),x2(k)/) ! ! if(size(w,2).gt.1) then ! ! call s1424(spldom%totdomweight,1,1,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat) ! ! else ! ! call s1424(spldom%totdomweight,0,0,point,cellk%leftknot(1),cellk%leftknot(2),wtmp,sstat) ! ! end if ! ! if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splinedomweight" ! ! if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splinedomweight" ! ! w(k,:)=wtmp(1:size(w,2)) ! ! ! CALL getgrad(spldom%totdomweightspl, x1(k), x2(k), wtmp(1), wtmp(2), wtmp(3)) ! w(k,:)=wtmp(1:size(w,2)) ! END SELECT !end DO if (size(w,2).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 !Do k=1,size(x2,1) ! SELECT CASE (spldom%cellk(i(k),j(k))%spltotkind) ! CASE(-1) ! w(k,:)=0 ! w(k,0)=-1 ! CYCLE ! CASE(1) ! w(k,:)=0 ! w(k,0)=1 ! CYCLE ! CASE DEFAULT ! end select !end do 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,sstat,sizew Integer,allocatable::i(:),j(:) !type(cellkind):: cellk allocate(gtmp(size(x2,1),0:size(g,2)-1)) allocate(i(size(x2,1)),j(size(x2,1))) call getindex(x1, x2, spldom, i, j) if(present(w)) then gtmp=w else !Do k=1,size(x2,1) ! ! cellk=spldom%cellk(i(k),j(k)) ! call s1424(spldom%Dirdomweight,1,1,(/x1(k),x2(k)/),cellk%leftknot(3),cellk%leftknot(4),gtmp(k,:),sstat) ! if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1424 for splineweight" ! if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1424 for splineweight" !end DO CALL getgrad(spldom%Dirdomweightspl, x1, x2, gtmp(1,:), gtmp(2,:), gtmp(3,:)) end if Do k=1,size(x2,1) !cellk=spldom%cellk(i(k),j(k)) if(spldom%cellk(i(k),j(k))%spldirkind.eq.0)then if(gtmp(k,0) .ge. 0) then if(size(g,2) .gt. 1) then g(k,1:2)=-gtmp(k,1:2)*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val end if g(k,0)=(1-gtmp(k,0))*spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val else g(k,0)=spldom%boundaries(spldom%cellk(i(k),j(k))%linkedboundaries(1))%Dirichlet_val if(size(g,2).gt. 1) then g(k,1:2)=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 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) 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 curvepos(3:4)=curvepos(3:4)/sqrt(curvepos(3)**2+curvepos(4)**2) ! 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 ) 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 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, sstat, 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 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 if(numintcu.ge.1) then posres=0.5*(intcurve(1)%epar1(1)+intcurve(1)%epar1(2)) end if 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,sstat, 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) 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) ! 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) end do x2(nbeval2-2)=(spldom%x2(spldom%nb2-1)+spldom%x2(spldom%nb2))/2.0_db x2(nbeval2-1)=spldom%x2(spldom%nb2) 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 !!$OMP PARALLEL DO private(i,cellk) !do i=0,spldom%nb1 ! j=spldom%nb2 ! cellk=spldom%cellk(min(i,spldom%nb1-1),min(j,spldom%nb2-1)) ! If(abs(cellk%spldirkind) .eq. 1) Then ! wpredir(1,i,j)=cellk%spldirkind ! else ! call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), spldom%x1(i),spldom%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)), spldom%x1(i),spldom%x2(j), wpretot(:,i,j),spldom%dist_extent) ! end IF !end do !!$OMP END PARALLEL DO !call s1537(reshape(wpredir(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%Dirdomweight,sstat) !if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight" !if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight" CALL set_splcoef((/3,3/),x1,x2,spldom%Dirdomweightspl) call get_dim(spldom%Dirdomweightspl,dims) 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) !call s1537(reshape(wpretot(1,:,:),(/(spldom%nb1+1)*(spldom%nb2+1)/)),spldom%nb1+1,spldom%nb2+1,1,spldom%x1,spldom%x2,0,0,0,0,4,4,1,1,spldom%totdomweight,sstat) !if (sstat > 0 ) WRITE(*,*) "Warning ",sstat," in distance calculation s1537 for splineweight" !if (sstat < 0 ) WRITE(*,*) "Error ",sstat," in distance calculation s1537 for splineweight" 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, ip, ii, jj, nidbas(2) DOUBLE PRECISION :: temp0(SIZE(xp),sp%sp2%order), temp1(SIZE(xp),sp%sp2%order) DOUBLE PRECISION, ALLOCATABLE, SAVE :: funx(:,:), funy(:,:) DOUBLE PRECISION, ALLOCATABLE, SAVE :: ftemp0(:), ftemp1(:) 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