diff --git a/src/geometry_mod.f90 b/src/geometry_mod.f90 index 24bedc7..e99d965 100644 --- a/src/geometry_mod.f90 +++ b/src/geometry_mod.f90 @@ -1,1436 +1,1436 @@ !------------------------------------------------------------------------------ ! 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 b-splines and !> can load the definition of the geometry from input files !------------------------------------------------------------------------------ MODULE geometry USE constants USE bsplines USE mumps_bsplines USE splinebound IMPLICIT NONE type innerspline Integer, Allocatable:: k(:) ! Index in reduced set real(kind=db), Allocatable:: weight(:) ! geomtric weight at relevant cell end type innerspline Real(kind=db),save:: z_0=0, r_0=0, invr_r=0, invr_z=0, r_a=0, r_b=0, z_r=1, r_r=1, alpha=0 Real(kind=db),save:: r_bLeft=0, r_bRight=0 Real(kind=db), save:: Phidown=0,Phiup=0 Real(kind=db), save:: L_r=0.16, L_z=0.24 Logical, save:: nlweb=.false. Integer, save::Interior=-1 Integer, save:: above1=1, above2=-1 Integer, save :: walltype = 0 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 Real(kind=db), save, allocatable:: gridwgeom(:,:) ! Stores the geometric 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 type(mumps_mat):: etildet ! Transpose of Matrix of extendend web splines integer,save :: nbreducedspline ! Number of splines in the reduced set type(spline_domain):: the_domain PROCEDURE(geom_eval), POINTER:: total_weight => NULL() PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() PUBLIC:: geom_weight ABSTRACT INTERFACE SUBROUTINE gtilde_eval(z,r,g) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(:,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), INTENT(OUT), OPTIONAL:: wupper(:,0:) END SUBROUTINE END INTERFACE INTERFACE geom_weight MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2 END INTERFACE geom_weight NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight contains SUBROUTINE read_geom(Fileid, rnorm, splrz, Potinn, Potout) Use mpi Use bsplines use basic, ONLY: phinorm type(spline2d):: splrz Real(kind=db):: rnorm, Potinn, Potout Integer:: Fileid, mpirank, ierr, istat character(len=1000) :: line !type(spline_boundary):: border0 !Real(kind=db):: cpoints(2,3),t(1), distance, pos !Real(kind=db), ALLOCATABLE::result(:,:,:) 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 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 defined 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 L_r=L_r/rnorm L_z=L_z/rnorm SELECT CASE (walltype) CASE(-3) ! test case total_gtilde=>gtest total_weight=>geom_w2 CASE(-2) ! test case with splinebound total_gtilde=>gtest total_weight=>geom_spline call read_splinebound(Fileid,the_domain, splrz, rnorm, Phinorm, Potinn,Potout) CASE(-1) ! test case total_gtilde=>gtest total_weight=>geom_weightstd CASE (2) ! coaxial cylinder and top ellipse with cylinder extensions total_gtilde=>gUpDown total_weight=>geom_w2 CASE (3) ! Two ellipses with "parallel" tangents with cylinders total_gtilde=>gUpDown total_weight=>geom_w3 CASE (4) ! Two ellipses with same radii with cylinders total_gtilde=>gUpDown total_weight=>geom_w4 CASE (5) ! Two ellipses with same radii with cylinders and total_gtilde=>gUpDown total_weight=>geom_w5 CASE (6) ! circular coaxial tilted ellipse right Dirichlet total_gtilde=>gUpDown total_weight=>geom_w6 CASE (7) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown total_weight=>geom_w7 CASE (8) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown total_weight=>geom_w8 CASE (9) - ! circular coaxial tilted ellipse right and left Dirichlet + ! Geometry defined as a spline curve total_gtilde=>gspline total_weight=>geom_spline call read_splinebound(Fileid,the_domain, splrz, rnorm, phinorm, Potinn,Potout) CASE DEFAULT ! Ellipse as in gt170 standard weight and straight coaxial configs total_gtilde=>gstd total_weight=>geom_weightstd END SELECT end subroutine read_geom SUBROUTINE init_geom(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 calssification for each cell Allocate(gridcelltype(nrz(1),nrz(2))) Call classifycell(zgrid, rgrid, gridcelltype) if (nlweb) then Allocate(bsplinetype(nrank(1)*nrank(2))) Call classifyspline(spl2, gridcelltype, bsplinetype) Call buildetilde(spl2,bsplinetype,gridcelltype,etilde) end if ALLOCATE(gridwgeom((nrz(1)+1)*(nrz(2)+1),0:2)) gridwgeom=0 ALLOCATE(gtilde((nrz(1)+1)*(nrz(2)+1),0:2)) gtilde=0 CALL total_gtilde(vec1, vec2, gtilde) Call geom_weight (vec1, vec2, gridwgeom) end Subroutine init_geom Subroutine geom_diag(File_handle, str, rnorm) Use mpi Use futils 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_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", L_r*RNORM) Call attach(File_handle, trim(grpname), "L_z", L_z*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',gridwgeom) END IF End subroutine geom_diag Subroutine buildetilde(spl2, bsplinetype, celltype, etilde) USE mumps_bsplines type(spline2d):: spl2 Integer:: bsplinetype(:) Integer:: celltype(:,:) type(mumps_mat):: etilde 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(:) !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(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(hmesh(nrank(1)*nrank(2))) allocate(d(size(bsplinetype))) Ibsplinetype=.False. 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 End do ! allocate etilde and its transpose call init(nrank(1)*nrank(2),nbreducedspline,etilde) call init(nrank(1)*nrank(2),nbreducedspline,etildet) k=1 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 !Ibsplinetype(i)=icellz-norder(1) .gt. 0 .and. icellr-norder(2) .gt. 0 .and. icellz .le. nrz(1) .and. icellr .le. nrz(2) 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 !!$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=(rgridmesh(j)-rgridmesh)**2+(zgridmesh(j)-zgridmesh)**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(indexdistance .gt. ((norder(1)+1)**2+(norder(2)+1)**2))then Write(*,*) 'Error on system conditioning, the number of radial or axial points is too low!' Write(*,*) 'Distance found: ', sqrt(indexdistance) !stop end if 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))) 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*(jcellr-icellr-m)/(n-m)*(jcellz-icellz-l)/(i-l) 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)) 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 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.100) !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)/2.lt.1e-12) 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 combine an ellipse function and an horizontal wall ! !--------------------------------------------------------------------------- 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 total_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(1:1,0:size(w,1)-1) ztmp=z rtmp=r call total_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 total_weight(z,r,w) ! End SUBROUTINE geom_weight2 SUBROUTINE geom_spline(z,r,w,wupper) Use splinebound, ONLY: spline_w Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) call spline_w(the_domain,z,r,w) End SUBROUTINE geom_spline SUBROUTINE geom_weightstd(z,r,w,wupper) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) Real(kind=db):: squareroot(size(w,1)), denom(size(w,1)) 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(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 geom_w2(z,r,w,wupper) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) call cyllweight(z,r,walltmp, r_a, above1) call ellipsewall(z,r,elliptmp,r_b,above2,r_0, z_0, invr_z, invr_r, Interior) if( present(wupper) ) then wupper=elliptmp w=walltmp else call Combine(elliptmp, walltmp, w, -1) end if End SUBROUTINE geom_w2 SUBROUTINE geom_w3(z,r,w,wupper) ! Defines the geometric weight for two facing ellipses of parallel boundaries with line prolongations !-----\__/------ ! !---\______/---- Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: belowtmp(size(w,1),0:size(w,2)-1), abovetmp(size(w,1),0:size(w,2)-1) ! Weight functions necessary for calculation of top and bottom boundaries ! Ellipse and bottom cylinder combination call ellipsewall(z,r,belowtmp,r_a,1,r_b, z_0, 1/(z_r+r_b-r_a), 1/(z_r+r_b-r_a), 1) ! Ellipse and top cylinder combination call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1) if( present(wupper) ) then wupper=abovetmp w=belowtmp else call Combine(belowtmp, abovetmp, w, -1) end if End SUBROUTINE geom_w3 SUBROUTINE geom_w4(z,r,w,wupper) ! Defines the geometric weight for two facing ellipses of same dimensions with line prolongations !-----\__/------ ! !-----\__/------ Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: belowtmp(size(w,1),0:size(w,2)-1), abovetmp(size(w,1),0:size(w,2)-1) ! Weight functions necessary for calculation of g ! Ellipse and bottom cylinder combination call ellipsewall(z,r,belowtmp,r_a,1,r_a, z_0, invr_z, invr_r, 1) ! Ellipse and top cylinder combination call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1) if( present(wupper) ) then wupper=abovetmp w=belowtmp else call Combine(belowtmp, abovetmp, w, -1) end if End SUBROUTINE geom_w4 SUBROUTINE geom_w5(z,r,w, wupper) ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w5 SUBROUTINE geom_w6(z,r,w, wupper) use basic, ONLY: zgrid, nz ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! tilted ellipse call discweight(z,r, w1, zgrid(nz),-1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w2 w=w1 else ! gives total weight call Combine(w1, w2, w, -1) end if End SUBROUTINE geom_w6 SUBROUTINE geom_w7(z,r,w, wupper) use basic, ONLY: zgrid, nz ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region ! and left and right dirichlet boundaries Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! Right Dirichlet call discweight(z,r, w1, zgrid(nz),-1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! Left Dirichlet call discweight(z,r, w1, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w7 SUBROUTINE geom_w8(z,r,w, wupper) use basic, ONLY: zgrid, nz ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region ! and left and right dirichlet boundaries Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! Left Dirichlet call discweight(z,r, w1, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, -1) ! Right Dirichlet call discweight(z,r, w1, zgrid(nz),-1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w8 Subroutine Boundary_limits(pts) Real(kind=db), Allocatable, INTENT(out) :: pts(:,:) If(allocated(pts))Deallocate(pts) allocate(pts(4,2)) pts(1,1:2)=(/z_0-invr_z,r_0/) pts(2,1:2)=(/z_0,r_0+invr_r/) pts(3,1:2)=(/z_0+invr_z,r_0/) pts(4,1:2)=(/z_0,r_0-invr_r/) End subroutine Boundary_limits SUBROUTINE classification(z, r, ctype) ! check if cell is fully inside, outside or on the boundary of the domain real(kind=db), INTENT(IN):: r(2), z(2) INTEGER, INTENT(OUT):: ctype Real(kind=db)::weights(5,1:1) 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 geom_weight(zeval,reval,weights) ctype=int(sign(1.0_db,weights(5,1))) If(weights(1,1)*weights(2,1) .le. 0 ) then ctype=0 return End If If(weights(1,1)*weights(3,1) .le. 0 ) then ctype=0 return End If If(weights(2,1)*weights(4,1) .le. 0 ) then ctype=0 return End If If(weights(3,1)*weights(4,1) .le. 0 ) then ctype=0 return End If end subroutine ! ################################################################## Subroutine calc_gauss(spl2, ngauss, i, j, xgauss, wgauss, gausssize, celltype) 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(:) Integer:: i,j, ngauss(2) Real(kind=db),Allocatable:: pts(:,:), 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, w, 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 !Call Boundary_Limits(pts) !Do k=1,size(pts,1) ! hasmaxpoint=hasmaxpoint .or. (pts(k,1) .ge. zgrid(i) .and. pts(k,1) .le. zgrid(i+1) & ! & .and. pts(k,2) .ge. rgrid(j) .and. pts(k,2) .le. rgrid(j+1)) !End do If(.not. hasmaxpoint) Then directiondown=1 directionup=1 Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown) Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup) select case ( directionup+directiondown) Case (0) ! The intersections are only on the left and right cell boundaries nbzpoints=2 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),zgrid(i+1)/) Case(1) nbzpoints=3 Allocate(zpoints(nbzpoints)) if(directiondown.eq.1) zpoints=(/zgrid(i), xptdown,zgrid(i+1)/) if(directionup .eq. 1) zpoints=(/zgrid(i), xptup,zgrid(i+1)/) Case(2) nbzpoints=3 Allocate(zpoints(nbzpoints)) call geom_weight(zgrid(i),rgrid(j),w) If(w.ge.0) zpoints=(/zgrid(1),min(xptdown,xptup),max(xptdown,xptup) /) If(w.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /) 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 geom_weight(zg(1),rgrid(j),w) rmin=rgrid(j) if (w .lt. 0) rmin = rgrid(j+1) Do k=1,ngauss(1) direction=2 Call Find_crosspoint((/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 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 End If 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) Real(kind=db):: x(2), y Real(kind=db):: xptm, xpt, temp Real(kind=db):: w, wold Integer, Intent(INOUT):: direction Integer:: i xptm=x(1) xpt=x(2) i=0 if(direction .eq.1) Call geom_weight(xptm,y,wold) if(direction .eq.2) Call geom_weight(y,xptm,wold) if(direction .eq.1) Call geom_weight(xpt,y,w) if(direction .eq.2) Call geom_weight(y,xpt,w) if(w*wold.gt.0) then direction=0 return End If Do while( i .le.100 .and. w*wold .lt. 0) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold if(direction .eq.1) Call geom_weight(xpt,y,w) if(direction .eq.2) Call geom_weight(y,xpt,w) xptm=xpt-w*(xpt-xptm)/(w-wold) wold=w temp=xptm xptm=xpt xpt=temp i=i+1 End do if(xpt .ge. x(2) .or. xpt .le. x(1)) direction=0 End Subroutine SUBROUTINE cyllweight(z,r,w, r_lim, above) Real(kind=db):: z(:),r(:),w(:,0:), r_lim Integer:: above w(:,0)=above*(r-r_lim) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=0 ! z derivative of w w(:,2)=above ! r derivative of w End If End subroutine cyllweight SUBROUTINE discweight(z,r,w, z_lim, right) Real(kind=db):: z(:),r(:),w(:,0:), z_lim Integer:: right w(:,0)=right*(z-z_lim) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=right ! z derivative of w w(:,2)=0 ! r derivative of w End If End subroutine discweight SUBROUTINE tiltedwall(z,r,w,r0,z0,alpha,left) Real(kind=db):: z(:),r(:),w(:,0:) Integer:: left Real(kind=db):: slope, r0, z0, alpha slope=tan(alpha) w(:,0)=(r-r0-slope*(z-z0)) If(size(w,2) .gt. 1) then ! first derivative w(:,1)= -slope ! z derivative of w w(:,2)= 1 ! r derivative of w End If w=left*w end SUBROUTINE SUBROUTINE ellipsewall(z,r,w,r_lim,above,r_0, z_0, invr_z, invr_r, Interior) Real(kind=db):: z(:),r(:),w(:,0:) Integer:: Interior, above Real(kind=db):: r_lim, r_0,z_0,invr_z,invr_r Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) call cyllweight(z,r,walltmp, r_lim, above) call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior) call combine( walltmp, elliptmp, w, Interior) END SUBROUTINE ellipsewall Subroutine ellipseweight(z,r,w, r_0, z_0, invr_z, invr_r, Interior) Real(kind=db):: z(:),r(:),w(:,0:) Real(kind=db):: r_0,z_0,invr_z,invr_r Integer:: Interior Real(kind=db):: D(size(r,1)) D=sqrt((r-r_0)**2*invr_r**2+(z-z_0)**2*invr_z**2) w(:,0)=Interior*(1-D) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=-(z-z_0)*invr_z**2/D*Interior ! z derivative of w w(:,2)=-(r-r_0)*invr_r**2/D*Interior ! r derivative of w End If End subroutine ellipseweight Subroutine tiltedellipse(z,r,w, r_0, z_0, invr_z, invr_r, Interior, alpha) Real(kind=db):: z(:),r(:),w(:,0:) Real(kind=db):: r_0,z_0,invr_z,invr_r, alpha, cosa, sina Real(kind=db):: deltar(size(r,1)), deltaz(size(z,1)) Integer:: Interior Real(kind=db):: D(size(r,1)) cosa=cos(alpha) sina=sin(alpha) deltar=(r-r_0) deltaz=(z-z_0) D=sqrt(((deltaz*cosa+deltar*sina)*invr_z)**2+((deltaz*sina-deltar*cosa)*invr_r)**2) w(:,0)=1-D ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=(-cosa*(deltaz*cosa+deltar*sina)*invr_z**2-sina*(deltaz*sina-deltar*cosa)*invr_r**2)/D ! z derivative of w w(:,2)=(-sina*(deltaz*cosa+deltar*sina)*invr_z**2+cosa*(deltaz*sina-deltar*cosa)*invr_r**2)/D ! r derivative of w End If w=w*Interior End subroutine tiltedellipse SUBROUTINE combine(w1, w2, w, union) ! Combines two weights w1 and w2 into w ! If union is 1 this defines the union of the geometric domains ! If union is -1 this defines the intersection of the geometric domains Real(kind=db):: w1(:,0:), w2(:,0:), w(:,0:) Integer:: union ! if 1 defines the union ow weights, if -1 defines intersection Real(kind=db):: squareroot(size(w,1)), denom(size(w,1)) denom=w1(:,0)**2+w2(:,0)**2 squareroot=union*sqrt(denom) w(:,0)=w1(:,0)+w2(:,0)+squareroot ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=w1(:,1)+w2(:,1)+(w1(:,1)*w1(:,0)+w2(:,1)*w2(:,0))/squareroot ! z derivative of w w(:,2)=w1(:,2)+w2(:,2)+(w1(:,2)*w1(:,0)+w2(:,2)*w2(:,0))/squareroot ! r derivative of w End If END SUBROUTINE Subroutine gstd(z,r,gtilde) ! standard g function added to rhs of poisson solver Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-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,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 SUBROUTINE gspline(z,r,g) Use splinebound, ONLY: spline_g Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(:,0:) call spline_g(the_domain,z,r,g) End SUBROUTINE gspline SUBROUTINE gUpDown(z,r,gtilde) ! standard g function added to rhs of poisson solver Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1) Real(kind=db):: denom(size(r,1)) ! Weight functions necessary for calculation of g call total_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) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: wtmp(size(gtilde,1),0:size(gtilde,2)-1) ! call total_weight(z,r,wtmp) gtilde(:,0)=(sin(2*pi*z/L_z)*cos(2*pi*r/L_r) + 2) If(size(gtilde,2) .gt. 1) then ! first derivative gtilde(:,1)=(2*pi/L_z*cos(2*pi*z/L_z)*cos(2*pi*r/L_r))*(1-wtmp(:,0))-wtmp(:,1)*gtilde(:,0) gtilde(:,2)=(-2*pi/L_r*sin(2*pi*z/L_z)*sin(2*pi*r/L_r))*(1-wtmp(:,0))-wtmp(:,2)*gtilde(:,0) End If gtilde(:,0)=gtilde(:,0)*(1-wtmp(:,0)) End subroutine Subroutine ftest(z,r,f) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT)::f(:,0:) ! f(:,0)=2*pi*sin(2*pi*z/L_z)*( & & 2*pi*cos(2*pi*r/L_r)*(1/L_z**2+1/L_r**2) & & +1/(r*L_r)*sin(2*pi*r/L_r) ) End Subroutine ftest subroutine Reducematrix(full,reduced) use mumps_bsplines type(mumps_mat):: full, reduced, tempmat Real(kind=db), allocatable:: temparray(:), rsltarray(:) Integer:: fullrank, reducedrank Integer:: i,j call to_mat(full) fullrank=full%rank reducedrank=nbreducedspline call init(reducedrank,2,reduced) call init(fullrank,2,tempmat) ! Compute tempmat=(gallerkin* etildet)' !$OMP parallel private(j), firstprivate(temparray,rsltarray), default(shared) allocate(temparray(fullrank),rsltarray(fullrank)) !$OMP do do i=1,fullrank call getrow(etilde,i,temparray) rsltarray=vmx_mumps_mat_loc(full,temparray) DO j=1,fullrank CALL putele_sploc(tempmat%mat%row(i), j, rsltarray(j), .false.) END DO end do !$OMP end do !$OMP barrier !$OMP single call to_mat(tempmat) !$OMP end single ! Compute reduced=(etilde*tempmat)' !$OMP do do i=1,reducedrank call getrow(etilde,i,temparray) rsltarray=vmx_mumps_mat_loc(tempmat,temparray) DO j=1,reducedrank CALL putele_sploc(reduced%mat%row(i), j, rsltarray(j), .false.) END DO end do !$OMP end do deallocate(temparray,rsltarray) !$OMP end parallel end subroutine !!############################################################################################### ! subroutines imported from mumps_bsplines and sparse lib to have thread-private routines SUBROUTINE putele_sploc(arow, j, val, nlforce_zero) ! ! Put (overwrite) 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 LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero ! TYPE(elt), TARGET :: pre_root TYPE(elt), POINTER :: t, p LOGICAL :: rmnode ! pre_root%next => arow%row0 ! pre_root is linked to the head of the list. t => pre_root ! ! Remove node which has zero val or not? ! But never create new node with zero val ! rmnode = .TRUE. IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero ! DO WHILE( ASSOCIATED(t%next) ) p => t%next IF( p%index .EQ. j ) THEN IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val! t%next => p%next arow%nnz = arow%nnz-1 arow%row0 => pre_root%next ! In case the head is altered DEALLOCATE(p) ELSE p%val = val END IF RETURN END IF IF( p%index .GT. j ) EXIT t => t%next END DO ! ! Never create new node with zero val ! IF(ABS(val).GT.EPSILON(0.0d0)) THEN ALLOCATE(p) p = elt(j, val, t%next) t%next => p arow%nnz = arow%nnz+1 arow%row0 => pre_root%next END IF END SUBROUTINE putele_sploc ! SUBROUTINE getcol_loc(amat, j, arr) ! ! Get a column from sparse matrix ! USE mumps_bsplines TYPE(mumps_mat), INTENT(in) :: amat INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(out) :: arr(:) INTEGER :: i ! DO i=amat%istart,amat%iend CALL getele_loc(amat, i, j, arr(i)) END DO END SUBROUTINE getcol_loc SUBROUTINE getele_loc(mat, i, j, val) ! ! Get element (i,j) of sparse matrix ! USE mumps_bsplines TYPE(mumps_mat) :: mat INTEGER, INTENT(in) :: i, j DOUBLE PRECISION, INTENT(out) :: val INTEGER :: iget, jget INTEGER :: k, s, e ! iget = i jget = j IF(mat%nlsym) THEN IF( i.GT.j ) THEN ! Lower triangular part iget = j jget = i END IF END IF ! val = 0.0d0 ! Assume zero val if outside IF( iget.LT.mat%istart .OR. iget.GT.mat%iend ) RETURN ! IF(ASSOCIATED(mat%mat)) THEN CALL getele(mat%mat, iget, jget, val) ELSE s = mat%irow(iget) - mat%nnz_start + 1 e = mat%irow(iget+1) - mat%nnz_start k = isearch(mat%cols(s:e), jget) IF( k.GE.0 ) THEN val =mat%val(s+k) ELSE val = 0.0d0 ! Assume zero val if not found END IF END IF END SUBROUTINE getele_loc !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION vmx_mumps_mat_loc(mat, xarr) RESULT(yarr) ! ! Return product mat*x ! Use mumps_bsplines TYPE(mumps_mat) :: mat DOUBLE PRECISION, INTENT(in) :: xarr(:) DOUBLE PRECISION :: yarr(SIZE(xarr)) DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0 CHARACTER(len=6) :: matdescra INTEGER :: n ! n = mat%rank ! !#ifdef MKL IF(mat%nlsym) THEN matdescra = 'sun' ELSE matdescra = 'g' END IF ! CALL mkl_dcsrmv('N', n, n, alpha, matdescra, mat%val, & & mat%cols, mat%irow(1), mat%irow(2), xarr, & & beta, yarr) !#else ! yarr = 0.0d0 ! DO i=1,n ! DO j=mat%irow(i), mat%irow(i+1)-1 ! yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j)) ! END DO ! IF(mat%nlsym) THEN ! DO j=mat%irow(i)+1, mat%irow(i+1)-1 ! yarr(mat%cols(j)) = yarr(mat%cols(j)) & ! & + mat%val(j)*xarr(i) ! END DO ! END IF ! END DO ! Write(*,*) "linear mult") !#endif ! END FUNCTION vmx_mumps_mat_loc END MODULE geometry \ No newline at end of file diff --git a/src/splinebound_mod.f90 b/src/splinebound_mod.f90 index 331ae61..20a0783 100644 --- a/src/splinebound_mod.f90 +++ b/src/splinebound_mod.f90 @@ -1,550 +1,557 @@ 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 + Real(kind=db):: epsge=1.0e-8 !< geometric resolution used for calculating distances 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, istat, mpirank, ierr - real(kind=db):: dist_extent, Potinn, Potout + real(kind=db):: dist_extent, Potinn, Potout, epsge 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 + namelist /spldomain/ nbsplines, dist_extent, h5fname, dr, Lz, ra, rb, epsge 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) 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) + call setspline_boundary(spldom%boundaries(1), cpoints, degree, Potinn, epsge) 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 setspline_boundary(spldom%boundaries(2), cpoints, degree, Potout, epsge) 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):: Dval, epsge 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), "epsge", epsge) 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) + Call setspline_boundary(spldom%boundaries(i),reshape(points,(/dim,posdim(1)/dim /)), order-1, Dval/phinorm, epsge) 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) + subroutine setspline_boundary(b_curve, cpoints, degree, D_val, epsge) 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 + Real(kind=db):: D_val, epsge 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 + b_curve%epsge=epsge 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 + + epsco=0 + epsge=b_curve%epsge + 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