diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index 18781e4..331ae61 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,535 +1,550 @@ MODULE splinebound USE constants USE bsplines USE forSISL, only: newCurve, freeCurve, freeIntCurve, writeSISLcurve, writeSISLpoints Use forSISLdata IMPLICIT NONE TYPE spline_boundary ! all curves assume right handedness to set which side of the curve is inside or outside type(SISLCurve):: curve !type(SISLIntCurve):: intcurve Real(kind=db):: Dirichlet_val !< Value for the dirichlet boundary condition created by this boundary END TYPE spline_boundary type spline_domain integer:: nb_splines type(spline_boundary), allocatable:: boundaries(:) Real(kind=db):: dist_extent=0.1 !< distance used for the merging with the plateau function for the weight type(cellkind), ALLOCATABLE:: cellk(:,:) type(spline2d), pointer:: splrz => null() Integer:: nb1 Integer:: nb2 end type spline_domain type cellkind integer:: splkind=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)=-1 !< stores the spline curve indices in the spline_domain of the spline boundaries that are at a distance lower than dist_extent integer:: uniontype=0 !< defines if the two boundaries are linked by 1: union, -1: intersection or 0, no joint integer:: leftknot(2)=0 !< for each linkedboundary saves an initial guess for calculating the distance end type cellkind CONTAINS subroutine read_splinebound(Fileid, spldom, splrz, rnorm, Phinorm, Potinn,Potout) + use mpi Integer:: Fileid type(spline_domain):: spldom type(spline2d):: splrz real(kind=db):: rnorm, phinorm - Integer:: nbsplines + Integer:: nbsplines, istat, mpirank, ierr real(kind=db):: dist_extent, Potinn, Potout - Character(len=128):: h5fname="" + Character(len=128):: h5fname="", line integer:: nbpoints real(kind=db),allocatable:: cpoints(:,:) real(kind=db):: dz, dr, ctr, Lz, ra, rb integer:: degree=4, i namelist /spldomain/ nbsplines, dist_extent, h5fname, dr, Lz, ra, rb + CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) REWIND(fileid) - READ(fileid,spldomain) + 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) if (.not. trim(h5fname)=='' ) then call setspline_domain(spldom, splrz, dist_extent, 0) call splinebound_readh5domain(trim(h5fname),spldom, rnorm, phinorm) call classifycells(spldom) return end if nbpoints=600 allocate(cpoints(2,nbpoints)) dist_extent=dist_extent/rnorm!5*(splrz%sp2%knots(1)-splrz%sp2%knots(0)) nbsplines=2 call setspline_domain(spldom, splrz, dist_extent, nbsplines) dz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/(nbpoints-1) Lz=Lz/rnorm ctr=(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0))/2 do i=1,nbpoints cpoints(1,i)=splrz%sp1%knots(0)+(i-1)*dz end do cpoints(2,:)=ra/rnorm call setspline_boundary(spldom%boundaries(1), cpoints, degree, Potinn) do i=1,nbpoints cpoints(1,i)=splrz%sp1%knots(0)+(nbpoints-i)*dz end do dr=dr/rnorm cpoints(2,:)=rb/rnorm do i=1,nbpoints if(((cpoints(1,i)-ctr)/Lz)**2.le.1) then cpoints(2,i)=rb/rnorm+dr*sqrt(1-(cpoints(1,i)-ctr)**2/Lz**2) endif end do call setspline_boundary(spldom%boundaries(2), cpoints, degree, Potout) call classifycells(spldom) end subroutine Subroutine splinebound_diag(File_handle, str, spldom) Use mpi Use futils Use basic, ONLY: rnorm, phinorm Integer:: File_handle type(spline_domain):: spldom Character(len=*):: str 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(str),"/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%nb_splines) do i=1,spldom%nb_splines Write(grpname,'(a,a,i2.2)') trim(str),"/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 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:: order,kind, dim INTEGER:: posrank, posdim(1) real(kind=db):: Dval real(kind=db),allocatable:: points(:) call openf(filename, h5id,'r') call getatt(h5id, '/geometry_spl/','nbsplines', spldom%nb_splines) call getatt(h5id, '/geometry_spl/','dist_extent', spldom%dist_extent) ! 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%nb_splines)) ! Read each boundary curve individually do i=1,spldom%nb_splines 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 Call getatt(h5id, trim(grpname), "Dirichlet_val", Dval) Call getatt(h5id, trim(grpname), "order", order) Call getatt(h5id, trim(grpname), "dim", dim) CALL getdims(h5id, TRIM(grpname)//"/pos", posrank, posdim) allocate(points(posdim(1))) CALL getarr(h5id, TRIM(grpname)//"/pos", points) points=points/rnorm Call setspline_boundary(spldom%boundaries(i),reshape(points,(/dim,posdim(1)/dim /)), order-1, Dval/phinorm) deallocate(points) end do end subroutine splinebound_readh5domain 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)) !Prepare structures to host singular spline boundaries spldom%nb_splines=nb_splines if(spldom%nb_splines.gt. 0) allocate(spldom%boundaries(nb_splines)) spldom%dist_extent=dist_extent end subroutine setspline_domain subroutine setspline_boundary(b_curve, cpoints, degree, D_val) Use bsplines use mpi type(spline_boundary):: b_curve Real(kind=db):: cpoints(:,:) Real(REAL64),ALLOCATABLE:: knots(:),points(:) integer:: degree Integer:: order, ierr Real(kind=db):: D_val Integer:: nbpoints,i, dim, kind, copy 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 kind=1 ! polynomial b-spline curve copy = 1 ! copy input array Allocate(knots(1:nbpoints+order)) DO i=order,nbpoints+1 knots(i)=real(i-order)/real(nbpoints+1-order) END DO knots(1:order)=knots(order) knots(nbpoints+1:nbpoints+order)=knots(nbpoints+1) allocate(points(dim*nbpoints)) points=reshape(cpoints,(/dim*nbpoints/)) CALL newCurve(nbpoints, order, knots, points, kind, dim, copy, b_curve%curve) b_curve%Dirichlet_val=D_val end subroutine setspline_boundary !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the 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_w(spldom,x1,x2,w) type(spline_domain):: spldom Real(kind=db), INTENT(IN):: x2(:),x1(:) Real(kind=db), INTENT(OUT):: w(:,0:) Integer:: i,j,k type(cellkind):: cellk Do k=1,size(x2,1) call getindex(x1(k), x2(k), spldom, i, j) cellk=spldom%cellk(i,j) If(cellk%splkind .ne. 0) Then w(k,:)=0 w(k,0)=cellk%splkind CYCLE end IF call splineweight(spldom%boundaries(cellk%linkedboundaries(1)), x1(k), x2(k), w(k,:),spldom%dist_extent,leftknot=cellk%leftknot(1)) end DO End SUBROUTINE spline_w !--------------------------------------------------------------------------- !> @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,z,r,g) type(spline_domain):: spldom Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(:,0:) Real(kind=db):: denom,gtmp(0:2) Integer:: k,i,j type(cellkind):: cellk type(spline_boundary):: splb !Do k=1,size(r,1) ! call getindex(z(k), r(k), spldom, i, j) ! cellk=spldom%cellk(i,j) ! If(cellk%splkind .ne. 0) Then ! g(k,:)=0 ! g(k,0)=spldom%centerg ! CYCLE ! end IF ! splb=spldom%boundaries(cellk%linkedboundaries(1)) ! call splineweight(splb,z(k),r(k),gtmp,spldom%dist_extent,leftknot=cellk%leftknot(1)) ! if(g(k,0) .ge. 0) then ! denom=1/(1+(spldom%nb_splines-1)*gtmp(0)) ! if(size(g,2).gt. 1) then ! g(k,1:2)=gtmp(1:2)*(spldom%nb_splines)*(spldom%centerg-splb%Dirichlet_val)*denom**2 ! end if ! g(k,0)=(gtmp(0)*splb%invDirichlet_val+splb%Dirichlet_val)*denom ! else ! g(k,0)=splb%Dirichlet_val ! if(size(g,2).gt. 1) then ! g(k,1:2)=gtmp(1:2)*(spldom%nb_splines)*(spldom%centerg-splb%Dirichlet_val) ! end if ! end if !end DO Do k=1,size(r,1) call getindex(z(k), r(k), spldom, i, j) cellk=spldom%cellk(i,j) If(cellk%splkind .ne. 0) Then g(k,:)=0 g(k,0)=0 CYCLE end IF splb=spldom%boundaries(cellk%linkedboundaries(1)) call splineweight(splb,z(k),r(k),gtmp,spldom%dist_extent,leftknot=cellk%leftknot(1)) if(g(k,0) .ge. 0) then g(k,1:2)=-gtmp(1:2)*splb%Dirichlet_val g(k,0)=(1-gtmp(0))*splb%Dirichlet_val else g(k,0)=splb%Dirichlet_val if(size(g,2).gt. 1) then g(k,1:2)=-gtmp(1:2)*splb%Dirichlet_val end if end if end DO End SUBROUTINE spline_g !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Evaluates the b-spline curve and its derivatives at the spline parameter value t [0, 1] !> @param[in] b_curve spline_boundary containing the spline curve parameters !> @param[in] t (:) array of positions along the spline curve !> @param[out] evalpos(:,:,:) array of points for each position t first index correspond to the index in t, second index are the coordinates, third index defines the order of derivation by t !> @param[in] fder matrix of boundary interpolations g with first index corresponding to the position and second index to the derivative !--------------------------------------------------------------------------- subroutine eval(b_curve,t,evalpos,fder) Use forSISL, ONLY: s1227 type(spline_boundary):: b_curve Real(kind=db):: t(:) Real(kind=db),allocatable:: evalpos(:,:,:) integer,optional:: fder integer:: left,i,der, stat, dim, j REAL(real64),allocatable:: result(:) if( present(fder)) then der=fder else der=0 end if dim=b_curve%curve%idim allocate(result(dim*(der+1))) if (allocated(evalpos)) deallocate(evalpos) allocate(evalpos(dim,size(t,1),0:der)) do i=1,size(t,1) call s1227(b_curve%curve,der,t(i),left,result,stat) do j=0,der evalpos(:,i,der)=result(j*dim+1:(j+1)*dim) end do end do END subroutine !--------------------------------------------------------------------------- !> @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, leftknot) Use forSISL, ONLY: s1227,s1221 type(spline_boundary):: b_curve Real(kind=db)::r,z Real(kind=db):: weight(0:) Real(kind=db),OPTIONAL:: distance Integer,OPTIONAL:: leftknot integer:: stat, der, left,siz real(kind=db):: h, d, tpos real(kind=real64):: curvepos(2*b_curve%curve%idim) weight=0 der=1 if(present(leftknot)) left=leftknot call dist(b_curve,(/z,r/),d,tpos) ! position and derivative wrt r,z call s1221(b_curve%curve,der,tpos,left,curvepos,stat) if (stat > 0 ) WRITE(*,*) "Warning ",stat," in distance calculation s1227 for splineweight at ", z, r if (stat < 0 ) WRITE(*,*) "Error ",stat," in distance calculation s1227 for splineweight at ", z, r 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 if ((-(z-curvepos(1))*curvepos(4)+(r-curvepos(2))*curvepos(3)) .lt. 0 ) weight(0)=-weight(0) siz=size(weight,1) if (size(weight,1).gt.1 .and. 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)=-curvepos(4)/h !weight(2)=+curvepos(3)/h end if if(present(distance)) distance=d if(present(leftknot)) leftknot=left end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Calculates the closest distance between the point and the 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 type(spline_boundary):: b_curve Real(kind=db):: point(:) real(kind=db):: distance Real(kind=db),optional::pos REAL(real64):: posres, epsco, epsge, intpar integer:: numintpt, stat, numintcu epsco=1.9e-9 epsge=1.0e-8 numintpt=0 stat=0 call s1957(b_curve%curve,point,b_curve%curve%idim,epsco,epsge,posres,distance,stat) if (stat > 1 ) WRITE(*,*) "Warning ",stat," in distance calculation s1227 for splineweight at ", point if (stat < 0 ) WRITE(*,*) "Error ",stat," in distance calculation s1227 for splineweight at ", point if (present(pos)) pos=posres END subroutine SUBROUTINE classify(x1, x2, ctype, linkedspline, leftsave, spldom) real(kind=db), INTENT(IN):: x2(2), x1(2) INTEGER, INTENT(OUT):: ctype INTEGER, INTENT(OUT):: linkedspline INTEGER, INTENT(OUT):: leftsave type(spline_domain)::spldom Real(kind=db)::w(0:2) Real(kind=db):: zeval(4),reval(4) Real(kind=db):: distance, dmin, inside integer:: i,k, left dmin=1000*spldom%dist_extent linkedspline=0 zeval=(/ x1(1),x1(2),x1(1),x1(2) /) reval=(/ x2(1),x2(1),x2(2),x2(2) /) do i=1,spldom%nb_splines do k=1,4 call splineweight(spldom%boundaries(i),zeval(k),reval(k),w,spldom%dist_extent,distance,left) if(w(1) .ne. 0 .or. w(2) .ne. 0) then linkedspline=i ctype=0 leftsave=left return end if if(distance .lt. dmin) then dmin=distance inside=w(0) end if end do end do ctype=sign(1,int(inside)) end subroutine subroutine classifycells(spldom) type(spline_domain):: spldom integer:: i,j do i=0,spldom%nb1-1 do j=0,spldom%nb2-1 call classify(spldom%splrz%sp1%knots(i:i+1),spldom%splrz%sp2%knots(j:j+1),spldom%cellk(i,j)%splkind,spldom%cellk(i,j)%linkedboundaries(1), spldom%cellk(i,j)%leftknot(1),spldom) end do end do end subroutine subroutine getindex(x1,x2,spldom, i, j) type(spline_domain):: spldom real(kind=db):: x1, x2 integer:: i,j, nb1, nb2 call locintv(spldom%splrz%sp1,x1, i) call locintv(spldom%splrz%sp2,x2, j) end subroutine 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