Page MenuHomec4science
Files F100575968

splinebound_mod.f90
No OneTemporary
Actions

File Metadata

Created
Fri, Jan 31, 21:12

splinebound_mod.f90
View Options

!------------------------------------------------------------------------------
! EPFL/Swiss Plasma Center
!------------------------------------------------------------------------------
!
! MODULE: splinebound
!
!> @author
!> G. Le Bars - EPFL/SPC
!
!> Last modif.
!> 06/07 2023
!
! DESCRIPTION:
!> Module handling domain boundaries defined by NURBS curves
!------------------------------------------------------------------------------
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
real(kind=db),allocatable:: slims(:,:) !< Spline parameter limits for boundary i in the cell
real(kind=db),allocatable:: blen(:) !< curve length boundary i in the cell
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:: istat, mpirank, ierr
real(kind=db):: dist_extent ! distance in [m] over which the w goes from 1 to 0
Character(len=128):: h5fname="" ! name of the h5 file storing the geometry
Character(len=128):: line
real(kind=db) :: Dvals(30)=0
integer:: i
namelist /spldomain/ 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 precomputecells(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, 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))
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_db ! starting parameter for the knots vector
! initialize a new curve using SISL
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,x2,x1,w)
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(2,:), w(1,:))
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,x2,x1,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,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(2,:), w(1,:))
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,x2,x1,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,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
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(2,:), gtmp(1,:))
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 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
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-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)*abs((h-d))/h**3
weight(2)=+3*curvepos(3)*abs((h-d))/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, 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 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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> evaluate the cell limited by x1 and x2 and classify if this is an inner outer or boundary spline
!> in addition, for boundary cells, it stores the closest spline curve for the evaluation of the boundary and the
!> dirichlet weight
!> @param[in] x1 axial limits of the cell
!> @param[in] x2 radial limits of the cell
!> @param[in] cellk structure storing the parameters of this cell
!> @param[in] spldom structure defining the domain boundaries
!> @param[out] wpredir dirichletweight at the center of the cell
!> @param[out] wpretot totweight at the center of the cell
!---------------------------------------------------------------------------
SUBROUTINE classify(x1, x2, cellk, spldom, wpredir, wpretot)
Use forSISL, ONLY: s1227
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):: dmindir,dmintot, insidedir, insidetot, distance
integer:: i,k, left, sstatus
logical:: isboundarycell
real(kind=db):: curvepos1(2), curvepos2(2)
allocate(guess(spldom%nbsplines,4))
allocate(w(0:2,spldom%nbsplines,4))
w=0
cellk%spldirkind=0
guess=-1.0_db
dmindir=HUGE(spldom%dist_extent)
dmintot=HUGE(spldom%dist_extent)
cellk%linkedboundaries=0
allocate(cellk%blen(spldom%nbsplines),cellk%slims(2,spldom%nbsplines))
cellk%blen=0
cellk%slims=0
! we define the evaluation positions
!
! x(4)------x(3)
! | |
! | |
! x(1)------x(2)
reval=(/ x1(1),x1(2),x1(2),x1(1) /)
zeval=(/ x2(1),x2(1),x2(2),x2(2) /)
insidedir=1
insidetot=1
do i=1,spldom%nbsplines
isboundarycell=.false.
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((spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar) .and. distance .lt. dmindir) 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
cellk%linkedboundaries(1)=i
cellk%spldirkind=0
end if
dmindir=distance
! Otherwise we define the interior by the closest spline
insidedir=w(0,i,k)
end if
if(distance .lt. dmintot) 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
cellk%linkedboundaries(2)=i
cellk%spltotkind=0
end if
dmintot=distance
! Otherwise we define the interior by the closest spline
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(w(0,i,k).lt.0)then
insidetot=w(0,i,k)
dmintot=0
if(distance.lt.spldom%dist_extent)then
cellk%linkedboundaries(2) =i
else
cellk%linkedboundaries(2) =0
cellk%spltotkind=-1
end if
end if
! The negative boundaries take precedence
! this is important when they define what is outside of the simulation domain.
if((spldom%boundaries(i)%type .eq. bd_Dirichletconst .or. spldom%boundaries(i)%type .eq.bd_Dirichletvar).and. w(0,i,k).lt.0)then
dmindir=0
insidedir=w(0,i,k)
if(distance.lt.spldom%dist_extent)then
cellk%linkedboundaries(1) =i
else
cellk%linkedboundaries(1) =0
cellk%spldirkind=-1
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
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> evaluate all the cells of the FEM and classify if they are an inner outer or boundary spline
!> in addition, for boundary cells, it stores the closest spline curve for the evaluation of the boundary and the
!> dirichlet weight
!> This function also prepare the bivariate spline interpolant for quickly evaluating
!> the dirichletweight and the total weight
!> @param[in] spldom structure defining the domain boundaries
!---------------------------------------------------------------------------
subroutine precomputecells(spldom)
use forSISL, ONLY: s1537, s1424
use bsplines
type(spline_domain):: spldom
integer:: i,j, dims(2), nbeval1, nbeval2,k,l, order(2), evalorder(2)
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
order=(/3,3/)
evalorder=(/0,0/)
! 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) collapse(2)
do j=0,spldom%nb2-1
do i=0,spldom%nb1-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) collapse(2)
do j=0,nbeval2-1
do i=0,nbeval1-1
call locintv(spldom%splrz%sp1,x1(i),k)
call locintv(spldom%splrz%sp2,x2(j),l)
cellk=spldom%cellk(k,l)
If(abs(spldom%cellk(k,l)%spldirkind) .eq. 1) Then
wpredir(1,i,j)=spldom%cellk(k,l)%spldirkind
else
call splineweight(spldom%boundaries(spldom%cellk(k,l)%linkedboundaries(1)), x2(j),x1(i), wpredir(:,i,j),spldom%dist_extent)
end IF
If(abs(spldom%cellk(k,l)%spltotkind) .eq. 1) Then
wpretot(1,i,j)=spldom%cellk(k,l)%spltotkind
else
call splineweight(spldom%boundaries(spldom%cellk(k,l)%linkedboundaries(2)), x2(j),x1(i), 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(order,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, evalorder, 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(order,x1,x2,spldom%totdomweightspl)
call get_splcoef(spldom%totdomweightspl, wpretot(1,:,:), c)
CALL gridval(spldom%totdomweightspl, spldom%x1(1), spldom%x2(1), val, evalorder, c)
deallocate(c)
end subroutine
!---------------------------------------------------------------------------
!> @author
!> Guillaume Le Bars EPFL/SPC
!
! DESCRIPTION:
!>
!> @brief
!> evaluate the FEM cell index for the quick evaluation of the dirichletweight and the totalweight
!> using the bivariated b-spline interpolant
!> @param[in] x1 array of axial positions
!> @param[in] x2 array of radial positions
!> @param[in] spldom structure defining the domain boundaries
!> @param[out] i array of axial cell indices
!> @param[out] j array of radial cell indices
!---------------------------------------------------------------------------
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
!---------------------------------------------------------------------------
!> @author
!> adapted from spclibs
!
! DESCRIPTION:
!>
!> @brief
!> modified function from spclibs to evaluate bivariate b-splines and its derivatives at the positions xp,yp
!> @param[in] sp spline2d structure storing the bivariate b-spline interpolant data
!> @param[in] xp array of axial positions
!> @param[in] yp array of radial positions
!> @param[in] leftx array of axial cell indices for xp
!> @param[in] lefty array of radial cell indices for yp
!> @param[out] f00 array of evaluated function at xp, yp
!> @param[out] f10 array of axial derivative of the evaluated function at xp, yp
!> @param[out] f01 array of radial derivative of the evaluated function at xp, yp
!---------------------------------------------------------------------------
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(4:)
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(4:)
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

Event Timeline

c4science · Help