!------------------------------------------------------------------------------ ! EPFL/Swiss Plasma Center !------------------------------------------------------------------------------ ! ! MODULE: geometry ! !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> Module responsible for handling geometries with non constant radius using b-splines interpolation !> This module defines ways to comupte the weight function needed for weighted extended b-splines and !> can load the definition of the geometry from input files !> This module is based on the theory by K. Hollig and book "Finite element methods with b-splines" !> SIAM Frontiers in applied mathematics !------------------------------------------------------------------------------ MODULE geometry USE constants USE bsplines USE mumps_bsplines USE splinebound use weighttypes IMPLICIT NONE type innerspline Integer, Allocatable:: k(:) ! Index in reduced set real(kind=db), Allocatable:: weight(:) ! geomtric weight at relevant cell end type innerspline type test_params !< parameters defining the manufactured test solutionwhen negative weighttypes is used real(kind=db):: z0 real(kind=db):: r0 real(kind=db):: Lz real(kind=db):: Lr end type real(kind=db), save:: testkr=1, testkz=1 Logical, save:: nlweb=.true. !< use weighted extended b-splines and not only weighted b-splies Integer, save :: walltype = 0 !< type of geometric weight to use (see readgeom) Integer, save, Allocatable:: bsplinetype(:) !< Array containing the inner/outer type for each bspline Integer, save, Allocatable:: gridcelltype(:,:) !< Array containing the inner/outer type for each gridcell Integer, save, Allocatable:: linkedspline(:,:) !< Array containing the lowerleft linked spline in case of boundary spline Real(kind=db), save, allocatable:: gridwdir(:,:) !< Stores the Dirichlet geometric weight at the grid points Real(kind=db), save, allocatable:: gridwdom(:,:) !< Stores the domain weight at the grid points Real(kind=db), save, allocatable:: gtilde(:,:) ! Stores the extension to the domain of the boundary conditions Real(kind=db), save, allocatable:: fsolution(:,:) ! Stores the manufactured solution type(mumps_mat):: etilde !> Matrix of extendend web splines definition 4.9 p48 of Hollig's book type(mumps_mat):: etildet ! Transpose of Matrix of extendend web splines integer,save :: nbreducedspline ! Number of splines in the reduced set type(test_params), save:: test_pars PROCEDURE(geom_eval), POINTER:: dirichlet_weight => NULL()!< Function evaluating the weight for Dirichelt boundary conditions PROCEDURE(geomtot_eval), POINTER:: domain_weight => NULL() !< function giving the limits of the simulation domain PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() !< Computes the parameter gtilde used to impose dirichlet boundary conditions with phi=uh+gtilde PUBLIC:: geom_weight, dom_weight ABSTRACT INTERFACE SUBROUTINE gtilde_eval(z,r,g,w) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) END SUBROUTINE SUBROUTINE geom_eval(z,r,w,wupper) USE constants Real(kind=db), INTENT(IN) :: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), OPTIONAL :: wupper(0:,:) END SUBROUTINE SUBROUTINE geomtot_eval(z,r,w,idwall) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) INTEGER, optional, INTENT(OUT):: idwall(:) END SUBROUTINE END INTERFACE INTERFACE geom_weight MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2 END INTERFACE geom_weight INTERFACE dom_weight MODULE PROCEDURE dom_weight0, dom_weight1, dom_weight2, dom_weight3 END INTERFACE dom_weight NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, z_a, z_b ,walltype, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight, testkr, testkz contains !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Read the input parameters to initialize the geometry module from the standard input file !> @param[in] Fileid Text file id of the input file containing namelists !> @param[in] rnorm distance normalization constant !> @param[in] splrz bspline structure used by the FEM comming form bspline library !> @param[in] Potinn Normalized electric potential on the inner boundary !> @param[in] Potout Normalized electric potential on the outer boundary !--------------------------------------------------------------------------- SUBROUTINE read_geom(Fileid, rnorm, splrz, Potinn, Potout) use mpi Use bsplines use basic, ONLY: phinorm use weighttypes type(spline2d):: splrz Real(kind=db):: rnorm ! normalisation variable for distances Real(kind=db):: Potinn ! potential at inner electrode from basic Real(kind=db):: Potout ! potential at outer electrode from basic Integer:: Fileid, mpirank, ierr, istat character(len=1000) :: line CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) Rewind(Fileid) READ(Fileid, geomparams, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(Fileid) read(Fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, geomparams) !! Normalizations and initialization of geometric variables r_a=r_a/rnorm r_b=r_b/rnorm z_a=z_a/rnorm z_b=z_b/rnorm r_bLeft=r_bLeft/rnorm r_bRight=r_bRight/rnorm if(r_a .eq. 0 .and. r_b .eq.0) then !! in case no geom_params have been definedwe take the defaults from the grid limits r_a=splrz%sp2%knots(0) r_b=splrz%sp2%knots(splrz%sp2%nints) end if z_0=z_0/rnorm r_0=r_0/rnorm z_r=z_r/rnorm r_r=r_r/rnorm invr_r=1/r_r invr_z=1/z_r Phidown=Potinn Phiup=Potout SELECT CASE (abs(walltype)) CASE (2) ! coaxial cylinder and top ellipse with cylinder extensions total_gtilde=>gUpDown Dirichlet_weight=>geom_w2 domain_weight=>geom_rvaschtot CASE (3) ! Two ellipses with "parallel" tangents with cylinders total_gtilde=>gUpDown Dirichlet_weight=>geom_w3 domain_weight=>geom_rvaschtot CASE (4) ! Two ellipses with same radii with cylinders total_gtilde=>gUpDown Dirichlet_weight=>geom_w4 domain_weight=>geom_rvaschtot CASE (5) ! Two ellipses with same radii with cylinders and total_gtilde=>gUpDown Dirichlet_weight=>geom_w5 domain_weight=>geom_rvaschtot CASE (6) ! circular coaxial tilted ellipse right Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w6 domain_weight=>geom_rvaschtot CASE (7) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w7 domain_weight=>geom_rvaschtot CASE (8) ! circular coaxial tilted ellipse right and left Dirichlet total_gtilde=>gUpDown Dirichlet_weight=>geom_w8 domain_weight=>geom_rvaschtot CASE (9) ! Geometry defined as a spline curve total_gtilde=>gspline Dirichlet_weight=>geom_spline domain_weight=>geom_splinetot call read_splinebound(Fileid,the_domain, splrz, rnorm, phinorm) CASE (10) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w10 CASE (11) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w11 CASE (12) ! square section disc total_gtilde=>gUpDown domain_weight=>geom_rvaschtot Dirichlet_weight=>geom_w12 CASE DEFAULT ! Ellipse as in gt170 standard weight and straight coaxial configs total_gtilde=>gstd Dirichlet_weight=>geom_weightstd domain_weight=>geom_rvaschtot END SELECT ! If we are lauching a test case, we load the test_pars variable used for the manufactured solution if(walltype.lt.0) then test_pars%Lr=(splrz%sp2%knots(splrz%sp2%nints)-splrz%sp2%knots(0))/testkr test_pars%Lz=(splrz%sp1%knots(splrz%sp1%nints)-splrz%sp1%knots(0))/testkz test_pars%r0=0.5*(splrz%sp2%knots(splrz%sp2%nints)+splrz%sp2%knots(0)) test_pars%z0=0.5*(splrz%sp1%knots(splrz%sp1%nints)+splrz%sp1%knots(0)) total_gtilde=>gtest end if end subroutine read_geom !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Initialises the module and precomputes the matrix e used for extended splines !> Classify every grid cell (inner, outer, boundary) !> @param[in] splr2 bspline structure used by the FEM comming form bspline library !> @param[in] vec1 Axial Meshgrid array to precompute the weigths on the crid points !> @param[in] vec2 Radial Meshgrid array to precompute the weigths on the crid points !--------------------------------------------------------------------------- SUBROUTINE geom_init(spl2, vec1, vec2) type(spline2d):: spl2 real(kind=db):: vec1(:),vec2(:) Real(kind=db), Allocatable:: zgrid(:),rgrid(:) Integer:: nrank(2), nrz(2) Call get_dim(spl2, nrank, nrz) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) ! create a table of the classification for each cell Allocate(gridcelltype(nrz(1),nrz(2))) Call classifycells(zgrid, rgrid, gridcelltype) if (nlweb) then ! if we use extended splines, we need to build the e matrix linking inner and outer splines Allocate(bsplinetype(nrank(1)*nrank(2))) Call classifyspline(spl2, gridcelltype, bsplinetype) Call buildetilde(spl2,bsplinetype,gridcelltype) end if ! Precompute the domain and dirichlet weights at the grid/cell positions ALLOCATE(gridwdir(0:2,(nrz(1)+1)*(nrz(2)+1))) gridwdir=0 ALLOCATE(gridwdom(0:2,(nrz(1)+1)*(nrz(2)+1))) gridwdom=0 ALLOCATE(gtilde(0:2,(nrz(1)+1)*(nrz(2)+1))) gtilde=0 ALLOCATE(fsolution(0:2,(nrz(1)+1)*(nrz(2)+1))) fsolution=0 Call geom_weight (vec1, vec2, gridwdir) Call dom_weight (vec1, vec2, gridwdom) ! Precompute the gtilde at the grid cell position CALL total_gtilde(vec1, vec2, gtilde, gridwdir) ! Precompute the manufactured solution at the grid cell position if (walltype .lt. 0) then CALL manufsolution(vec1, vec2, fsolution) end if end Subroutine geom_init !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Save this module run parameters to the h5 result file in the correct group !> @param[in] File_handle h5 file id of the result file !> @param[in] parentgroup h5 parent group where to save the simulation parameters !> @param[in] rnorm distance normalization constant !--------------------------------------------------------------------------- Subroutine geom_diag(File_handle, parentgroup, rnorm) use mpi Use futils use weighttypes Integer:: File_handle Real(kind=db):: rnorm Character(len=*):: parentgroup 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(parentgroup),"/geometry" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "r_a", r_a*RNORM) Call attach(File_handle, trim(grpname), "r_b", r_b*RNORM) Call attach(File_handle, trim(grpname), "z_a", z_a*RNORM) Call attach(File_handle, trim(grpname), "z_b", z_b*RNORM) Call attach(File_handle, trim(grpname), "z_0", z_0*RNORM) Call attach(File_handle, trim(grpname), "r_0", r_0*RNORM) Call attach(File_handle, trim(grpname), "r_r", r_r*RNORM) Call attach(File_handle, trim(grpname), "z_r", z_r*RNORM) Call attach(File_handle, trim(grpname), "L_r", test_pars%Lr*RNORM) Call attach(File_handle, trim(grpname), "L_z", test_pars%Lz*RNORM) Call attach(File_handle, trim(grpname), "interior", interior) Call attach(File_handle, trim(grpname), "above1", above1) Call attach(File_handle, trim(grpname), "above2", above2) Call attach(File_handle, trim(grpname), "walltype", walltype) Call putarr(File_handle, trim(grpname)//'/geomweight',transpose(gridwdom)) Call putarr(File_handle, trim(grpname)//'/dirichletweight',transpose(gridwdir)) Call putarr(File_handle, trim(grpname)//'/gtilde',transpose(gtilde)) Call putarr(File_handle, trim(grpname)//'/ctype',gridcelltype) Call putarr(File_handle, trim(grpname)//'/linked_s',linkedspline) Call putarr(File_handle, trim(grpname)//'/bsplinetype',bsplinetype) Call putarr(File_handle, trim(grpname)//'/fsolution',transpose(fsolution)) END IF End subroutine geom_diag !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Construct the e matrix used to link inner and boundary b-splines !> @param[in] spl2 bi-variate b-spline structure !> @param[in] bsplinetype 1D array classifying the individual b-splines as inner, boundary, outer !> @param[in] celltype 2D array classifying the individual grid cells as inner, boundary, outer !--------------------------------------------------------------------------- Subroutine buildetilde(spl2, bsplinetype, celltype) USE mumps_bsplines Use basic, ONLY: rnorm, mpirank type(spline2d):: spl2 Integer:: bsplinetype(:) Integer:: celltype(1:,1:) Logical, Allocatable:: Ibsplinetype(:) Integer:: i,j,k, icellz, icellr, jcellz, jcellr, nrank(2), nrz(2), norder(2) real(kind=db), allocatable:: zgrid(:), rgrid(:) real(kind=db):: wgeomi, indexdistance, eij integer:: l,m, n, linkedi type(innerspline):: innersplinelist real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:) integer, allocatable:: igridmesh(:), jgridmesh(:) !if(allocated(etilde)) deallocate(etilde) nbreducedspline=count(bsplinetype .eq. 1) Call get_dim(spl2,nrank,nrz,norder) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) allocate(linkedspline(nrank(1),nrank(2))) allocate(innersplinelist%k(nrank(1)*nrank(2)),innersplinelist%weight(nrank(1)*nrank(2))) allocate(Ibsplinetype(nrank(1)*nrank(2))) allocate(rgridmesh(nrank(1)*nrank(2)), zgridmesh(nrank(1)*nrank(2))) allocate(igridmesh(nrank(1)*nrank(2)), jgridmesh(nrank(1)*nrank(2))) allocate(hmesh(nrank(1)*nrank(2))) allocate(d(size(bsplinetype))) Ibsplinetype=.False. ! Compute the center of the 2D splines for the Lagrange interpolation call calcsplinecenters(spl2%sp1,cz,hz) call calcsplinecenters(spl2%sp2,cr,hr) Do i=0,nrank(2)-1 zgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cz rgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cr(i+1) hmesh(i*nrank(1)+1:(i+1)*nrank(1))=hr(i+1)*hz igridmesh(i*nrank(1)+1:(i+1)*nrank(1))=(/ (j,j=0,nrank(1)-1)/) jgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=i End do ! allocate memory for etilde and its transpose call init(nrank(1)*nrank(2),nbreducedspline,etilde) call init(nrank(1)*nrank(2),nbreducedspline,etildet) k=1 ! Compute the terms eii for the inner b-splines Do i=1,nrank(1)*nrank(2) if(bsplinetype(i) .ne. 1) cycle ! span of this spline is almost completely outside D ! We consider only splines with ! one cell of bspline i completely in D icellz=mod(i-1,nrank(1))+1 icellr=(i-1)/(nrank(1))+1 wgeomi=1 outer: do l=max(1,icellz-norder(1)),min(nrz(1),icellz) do m=max(1,icellr-norder(2)),min(nrz(2),icellr) if (celltype(l,m).eq.1) then call geom_weight((zgrid(l)+zgrid(l+1))/2,(rgrid(m)+rgrid(m+1))/2,wgeomi) EXIT outer end if end do end do outer call putele(etilde,k,i,1/wgeomi) call putele(etildet,i,k,1/wgeomi) innersplinelist%k(i)=k innersplinelist%weight(i)=1/wgeomi k=k+1 linkedspline(icellz,icellr)=i Ibsplinetype(i)=icellz+norder(1) .le. nrank(1) .and. icellr+norder(2) .le. nrank(2) if(.not. Ibsplinetype(i)) cycle do m=0,norder(2) do l=0,norder(1) ! Check if all positive splines in this spline domain are inner splines Ibsplinetype(i)=Ibsplinetype(i) .and. (bsplinetype(i+l+m*nrank(1)) .eq. 1) end do end do end do ! Compute the terms eij for the outer b-splines !!$OMP PARALLEL DO DEFAULT(SHARED), PRIVATE(d,k,icellz,icellr,jcellr,jcellz,indexdistance,n,i,m,l,eij,linkedi) Do j=1,nrank(1)*nrank(2) ! find the closest b-spline fully in D for the web method if(bsplinetype(j) .ne. 0) cycle d=0 where(Ibsplinetype)! calculate distance between center of Interior splines and spline j d=(zgridmesh(j)-zgridmesh)**2+(rgridmesh(j)-rgridmesh)**2 end where k=minloc(d,1,MASK=Ibsplinetype) icellz=mod(k-1,nrank(1)) icellr=(k-1)/(nrank(1)) jcellz=mod(j-1,nrank(1)) jcellr=(j-1)/(nrank(1)) indexdistance=real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db) if(d(k) .gt. ((norder(1)+2)**2+(norder(2)+2)**2)*2 .and. mpirank .eq. 0)then Write(*,'(a)') 'Warning on system conditioning, the number of radial or axial points could be too low!' Write(*,'(a,1f6.2,a,2(1pe12.4))') 'Distance found: ', sqrt(indexdistance), ' at (z,r): ',zgridmesh(j)*rnorm, rgridmesh(j)*rnorm !stop end if ! Compute the Lagrange polynomia linking spline i and spline j linkedspline(jcellz+1,jcellr+1)=k do n=0,norder(2) do i=0,norder(1) eij=1 !eij=1 do m=0,norder(2) if(n.eq.m) cycle eij=eij*(rgridmesh(j)-rgridmesh(k+m*nrank(1)))/(rgridmesh(k+n*nrank(1))-rgridmesh(k+m*nrank(1))) !eij=eij*real((jcellr-icellr-m),db)/real((n-m),db) end do do l=0,norder(1) if( i .eq. l ) cycle eij=eij*(zgridmesh(j)-zgridmesh(k+l))/(zgridmesh(k+i)-zgridmesh(k+l)) !eij=eij*real((jcellz-icellz-l),db)/real((i-l),db) end do linkedi=innersplinelist%k(k+i+n*nrank(1)) ! equivalent to findloc, necessary for ifort 17 eij=eij*innersplinelist%weight(k+i+n*nrank(1)) ! add the polynomia to the etilde matrix call putele(etilde,linkedi,j,eij) call putele(etildet,j,linkedi,eij) end do end do end do !!$OMP End parallel do call to_mat(etilde) call to_mat(etildet) end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Routine to compute the center of the 2d b-splines used in Lagrange interpolation !> @param[in] spl 1D b-spline structure !> @param[out] ctrs 1D array of spline centers !> @param[out] heights 1D array of the maximum amplitude of each spline function !--------------------------------------------------------------------------- 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(spl%knots(i-order)) x2=spl%knots(i-1)-1e5*Epsilon(spl%knots(i-1)) left1=max(0,i-order) left2=min(nx-2,i-2) call basfun(x1,spl,fun1, left1+1) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold call basfun(x2,spl,fun2, left2+1) fun3=1 j=0 Do while( j .le. 300) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold x3=(x1+x2)/2 call locintv(spl, x3, left3) call basfun(x3,spl,fun3, left3+1) if( (x2-x1).lt.1e-13) exit if(abs(fun3(i-left3,1)).lt.1e-15) exit if(fun3(i-left3,1)*fun1(i-left1,1).le.0) then fun2=fun3 x2=x3 left2=left3 else fun1=fun3 x1=x3 left1=left3 end if j=j+1 End do ctrs(i)=x3 heights(i)=fun3(i-left3,0) end do end if end subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Classify each grid cell as inner (1), boundary (0), or outer (-1) cell and store it in the !> gridctype 2D array !> @param[in] zgrid array of axial limits of the grid-cells !> @param[in] rgrid array of radial limits of the grid-cells !> @param[out] gridctype 2D array classifying the individual grid cells as inner, boundary, outer !--------------------------------------------------------------------------- Subroutine classifycells(zgrid,rgrid,gridctype) Real(kind=db):: zgrid(1:),rgrid(1:) Integer:: gridctype(:,:) Integer::i,j gridctype=-1 !$OMP parallel do private(i,j) Do j=1,size(rgrid,1)-1 Do i=1,size(zgrid,1)-1 ! Determines the type inner/boundary/outer for each cell Call classifycell((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),& & gridctype(i,j)) End Do End do !$OMP end parallel do End subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> For each spline i determines if its support is inside(1), partially inside (0) or outside (-1) of the simulation domain !> @param[in] spl2 structure storing the bi-variate b-splines !> @param[in] gridctype 2D array classifying the individual grid cells as inner, boundary, outer !> @param[out] bsptype array of b-spline classification !--------------------------------------------------------------------------- subroutine classifyspline(spl2,gridctype,bsptype) type(spline2d):: spl2 Integer:: gridctype(:,:) Integer:: bsptype(:) Integer:: nrank(2), nrz(2), ndegree(2) Integer:: i, j, mu, imin, imax, jmin, jmax Integer, Allocatable:: splinespan(:,:) call get_dim(spl2, nrank, nrz, ndegree) ! by default, all splines have part of their support outside the domain D bsptype=0 Allocate(splinespan(ndegree(1)+1,ndegree(2)+1)) Do mu=1,nrank(1)*nrank(2) ! scan the spline space ! obtain the axial spline index i=mod(mu-1,nrank(1))+1 ! obtain the radial spline index j=(mu-1)/nrank(1)+1 ! by default all cells are outside for correct behavior in boundaries splinespan=-1 ! find the axial span of this spline in cell indices imin=max(1,i-ndegree(1)) imax=min(nrz(1),i) ! find the radial span of this spline in cell indices jmin=max(1,j-ndegree(2)) jmax=min(nrz(2),j) ! obtain the cell type on which the spline is defined splinespan(1:(imax-imin+1),1:(jmax-jmin+1))=gridctype(imin:imax,jmin:jmax) if ( ANY( splinespan==1 ) ) then ! if at least one cell is fully in the domain the spline is an inside spline bsptype(mu)=1 else if (.not. ANY( splinespan==0 )) then ! if all the cells are outside, this is an outside spline bsptype(mu)=-1 end if end do End subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> general placeholder function used in fields_mod to compute the weight for weighted_bsplines !> This functions call the correct subroutine pointed by dirichlet_weight ! !> @param[in] z axial position !> @param[in] r radial position !> @param[out] w geometric weight !--------------------------------------------------------------------------- SUBROUTINE geom_weight0(z,r,w) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w Real(kind=db):: rtmp(1:1), ztmp(1:1), wtmp(1:1,1:1) ztmp=z rtmp=r call Dirichlet_weight(ztmp,rtmp,wtmp) w=wtmp(1,1) End SUBROUTINE geom_weight0 SUBROUTINE geom_weight1(z,r,w) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w(0:) Real(kind=db):: rtmp(1:1), ztmp(1:1) Real(kind=db):: wtmp(0:size(w,1)-1,1:1) ztmp=z rtmp=r call Dirichlet_weight(ztmp,rtmp,wtmp) w=wtmp(:,1) End SUBROUTINE geom_weight1 SUBROUTINE geom_weight2(z,r,w) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) call Dirichlet_weight(z,r,w) ! End SUBROUTINE geom_weight2 !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> general placeholder function used in fields_mod to compute the domain weight for weighted_bsplines !> This functions call the correct subroutine pointed by domain_weight !> if the weight is negative return the id of the closest wall. ! !> @param[in] z axial position !> @param[in] r radial position !> @param[out] w geometric weight !> @param[out] idwall unique identifier of the boundary on which the particle is lost !--------------------------------------------------------------------------- SUBROUTINE dom_weight0(z,r,w,idwall) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w Real(kind=db):: rtmp(1:1), ztmp(1:1), wtmp(1:1,1:1) INTEGER:: idwalltmp(1:1) INTEGER, optional, INTENT(OUT):: idwall ztmp=z rtmp=r call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp) w=wtmp(1,1) if (present(idwall)) idwall=idwalltmp(1) End SUBROUTINE dom_weight0 SUBROUTINE dom_weight1(z,r,w,idwall) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w(0:) Real(kind=db):: rtmp(1:1), ztmp(1:1) Real(kind=db):: wtmp(0:size(w,1)-1,1:1) INTEGER:: idwalltmp(1:1) INTEGER, optional, INTENT(OUT):: idwall ztmp=z rtmp=r call domain_weight(ztmp,rtmp,wtmp,idwall=idwalltmp) w=wtmp(:,1) if (present(idwall)) idwall=idwalltmp(1) End SUBROUTINE dom_weight1 SUBROUTINE dom_weight2(z,r,w,idwall) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) INTEGER, optional, INTENT(OUT):: idwall(:) call domain_weight(z,r,w,idwall=idwall) ! End SUBROUTINE dom_weight2 SUBROUTINE dom_weight3(z,r,w,idwall) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:) Real(kind=db):: wtmp(1:1,size(w,1)) INTEGER, optional, INTENT(OUT):: idwall(:) call domain_weight(z,r,wtmp,idwall=idwall) !Write(*,*) 'wtmp, size', wtmp, size(wtmp) w=wtmp(1,:) ! End SUBROUTINE dom_weight3 SUBROUTINE geom_weightstd(z,r,w,wupper) ! return the geometric weight for a coaxial configuration ! or central conductor + ellipse if walltype ==1 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), OPTIONAL:: wupper(0:,:) Real(kind=db):: walltmp(0:size(w,1)-1,size(w,2)), elliptmp(0:size(w,1)-1,size(w,2)) Real(kind=db):: squareroot(size(w,2)), denom(size(w,2)) call cyllweight(z,r,walltmp, r_a, above1) SELECT CASE (walltype) CASE (0) call cyllweight(z,r,elliptmp, r_b, above2) CASE DEFAULT call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior) END SELECT if (present(wupper))then w=walltmp wupper=elliptmp return end if if(interior.eq.0 .and. above2.eq.0) then w=walltmp(:,:) else if (above1 .eq. 0) then w=elliptmp(:,:) else denom=walltmp(0,:)**2+elliptmp(0,:)**2 squareroot=sqrt(denom) w(0,:)=walltmp(0,:)+elliptmp(0,:)-squareroot ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(1,:)=walltmp(1,:)+elliptmp(1,:)-(elliptmp(1,:)*elliptmp(0,:)+walltmp(1,:)*walltmp(0,:))/& & squareroot ! z derivative of w w(2,:)=walltmp(2,:)+elliptmp(2,:)-(elliptmp(2,:)*elliptmp(0,:)+walltmp(2,:)*walltmp(0,:))/& & squareroot ! r derivative of w End If End if End SUBROUTINE geom_weightstd !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Classify each grid cell as inner (1), boundary (0), or outer (-1) cell and store it in the !> gridctype 2D array !> @param[in] zgrid array of axial limits of the grid-cells !> @param[in] rgrid array of radial limits of the grid-cells !> @param[out] ctype classification of the cell !--------------------------------------------------------------------------- SUBROUTINE classifycell(z, r, ctype) ! classify if cell is fully inside, outside or on the boundary of the domain ! by calculating the weight on each corner and the cell. ! It is assumed that the cells are sufficiently small such that there is no sharp edge entering the cell ! where an inner portion of only one cell edge is outside of the domain real(kind=db), INTENT(IN):: r(2), z(2) INTEGER, INTENT(OUT):: ctype Real(kind=db)::weights(1:1,5) Real(kind=db):: zeval(5),reval(5) zeval=(/ z(1),z(2),z(1),z(2), (z(2)+z(1))/2 /) reval=(/ r(1),r(1),r(2),r(2), (r(2)+r(1))/2 /) CAll dom_weight(zeval,reval,weights) ctype=int(sign(1.0_db,weights(1,5))) If(weights(1,1)*weights(1,2) .le. 0 ) then ctype=0 return End If If(weights(1,1)*weights(1,3) .le. 0 ) then ctype=0 return End If If(weights(1,2)*weights(1,4) .le. 0 ) then ctype=0 return End If If(weights(1,3)*weights(1,4) .le. 0 ) then ctype=0 return End If If(weights(1,3)*weights(1,2) .le. 0 ) then ctype=0 return End If If(weights(1,1)*weights(1,4) .le. 0 ) then ctype=0 return End If end subroutine ! ################################################################## !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> calculates the gauss quadrature integration points for the FEM method for any cell type !> takes care of boundary cells as well and limit the integration boundaries accordingly !> returns the gauss quadrature points and weights for cell (i,j) !> @param[in] spl2 structure storing the bi-variate b-spline !> @param[in] ngauss 2D array of number of gauss ponts in z and r !> @param[in] i axial cell index !> @param[in] j radial cell index !> @param[out] xgauss 2D array of evaluation points z=xgauss(:,1), r=xgauss(:,2) !> @param[out] xgauss 1D array of evaluation weights !> @param[out] gausssize total number of gauss points for this cell !> @param[out] celltype return the type of the cell inner (1), boundary (0), outer(-1) !--------------------------------------------------------------------------- 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, intent(out)::xgauss(:,:), wgauss(:) Integer:: i,j, ngauss(2) Real(kind=db),Allocatable:: zpoints(:) Real(kind=db):: zg(ngauss(1)),rg(ngauss(2)), wzg(ngauss(1)), wrg(ngauss(2)) Integer:: k, l, directiondown, directionup, nbzpoints, direction, ctype Logical:: hasmaxpoint Real(kind=db):: xptup, xptdown,wlu,wld, rmin, rmax, zmin, zmax 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(PRESENT(celltype)) celltype=ctype If (ctype .ge. 1) then ! we have a normal internal cell Allocate(xgauss(ngauss(1)*ngauss(2),2)) Allocate(wgauss(ngauss(1)*ngauss(2))) gausssize=ngauss(1)*ngauss(2) !Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO Else If(ctype.eq.0) then ! we have a boundary cell directiondown=1 directionup=1 ! We check if the boundary goes through the cell upper and lower limit Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown) Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup) call dom_weight(zgrid(i),rgrid(j),wld) call dom_weight(zgrid(i),rgrid(j+1),wlu) select case ( directionup+directiondown) Case (0) ! The intersections are only on the left and right cell boundaries ! or the upper and lower limits are full boundaries nbzpoints=2 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),zgrid(i+1)/) Case(1) if(directiondown.eq.1)then if( wlu .gt. 0) then ! the lower left corner is inside nbzpoints=3 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i), xptdown,zgrid(i+1)/) else nbzpoints=2 Allocate(zpoints(nbzpoints)) if(wld.gt.0)then ! the upper left corner is inside zpoints=(/zgrid(i),xptdown/) else zpoints=(/xptdown,zgrid(i+1)/) end if end if else if(wld .gt. 0) then ! the lower left corner is inside nbzpoints=3 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i), xptup,zgrid(i+1)/) else nbzpoints=2 Allocate(zpoints(nbzpoints)) if(wlu.gt.0)then ! the upper left corner is inside zpoints=(/zgrid(i),xptup/) else zpoints=(/xptup,zgrid(i+1)/) end if end if end if Case(2) nbzpoints=2 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i), zgrid(i+1)/) !If(wld.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /) !else ! nbzpoints=2 ! Allocate(zpoints(nbzpoints)) ! If(wld.ge.0) zpoints=(/zgrid(i),min(xptdown,xptup) /) ! If(wld.lt.0) zpoints=(/max(xptdown,xptup), zgrid(i+1) /) !end if Allocate(xgauss(ngauss(1)*ngauss(2),2)) Allocate(wgauss(ngauss(1)*ngauss(2))) gausssize=ngauss(1)*ngauss(2) ! Compute gauss points !CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) !Call gauleg(zpoints(l),zpoints(l+1),zg,wzg,ngauss(1)) call gauleg(rgrid(j),rgrid(j+1),rg,wrg,ngauss(2)) Do k=1,ngauss(2) ! We test if the lower or upper side is in the domain call dom_weight(zgrid(i),rg(k),wld) zmin=zgrid(i) if (wld .lt. 0) zmin = zgrid(i+1) direction=1 Call Find_crosspointdico((/zgrid(i),zgrid(i+1)/),rg(k),zmax,direction) ! We compute the radial limits at each z position Call gauleg(min(zmin,zmax),max(zmin,zmax),zg,wzg,ngauss(1)) ! We obtain the gauss w and pos for these boundaries if(direction .eq. 0.and. wld .lt. 0) then wzg=0 end if xgauss(k : ngauss(2)*ngauss(1) : ngauss(1),1)= zg xgauss(k : ngauss(2)*ngauss(1) : ngauss(1),2)= rg(k) wgauss(k : ngauss(2)*ngauss(1) : ngauss(1)) = wzg*wrg(k) End do return 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)) Do k=1,ngauss(1) ! We test if the lower or upper side is in the domain call dom_weight(zg(k),rgrid(j),wld) rmin=rgrid(j) if (wld .lt. 0) rmin = rgrid(j+1) direction=2 Call Find_crosspointdico((/rgrid(j),rgrid(j+1)/),zg(k),rmax,direction) ! We compute the radial limits at each z position Call gauleg(min(rmin,rmax),max(rmin,rmax),rg,wrg,ngauss(2)) ! We obtain the gauss w and pos for these boundaries if(direction .eq. 0.and. wld .lt. 0) then wrg=0 end if xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),1) = zg(k) xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),2) = rg wgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1)) = wrg*wzg(k) End do End Do Else gausssize=0 Allocate(xgauss(1:1,2)) Allocate(wgauss(1:1)) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration xgauss(1,1)=(zgrid(i)+zgrid(i+1))*0.5 xgauss(1,2)=(rgrid(j)+rgrid(j+1))*0.5 wgauss(1)=0 End If End Subroutine calc_gauss !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> computes quickly if a particle is inside the domain D and returns the id of the closest boundary !> if the particle is outside the domain !> @param[in] z axial particle position !> @param[in] r radial particle position !> @param[in] i axial cell index !> @param[in] j radial cell index !> @param[out] idwall id of the closest boundary if the particle is outside the domain !> @param[out] inside .true. if the particle is inside D !--------------------------------------------------------------------------- subroutine is_insidegeom(z,r,i,j,idwall,inside) REAL(kind=db), intent(in):: z,r INTEGER,INTENT(in):: i,j INTEGER:: idwall logical:: inside Real(kind=db):: weight if(i.gt.size(gridcelltype,1)-1.or. i.lt.0 .or. j.gt.size(gridcelltype,2)-1.or. j.lt.0 )then inside=.false. idwall=0 return end if select case(gridcelltype(i+1,j+1)) case(-1) inside=.false. idwall=0 case(1) inside=.true. idwall=0 case(0) call dom_weight(z, r, weight, idwall=idwall) inside=weight.gt.0 end select end subroutine is_insidegeom !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> calculates the boundary limit between x(1) and x(2) using Newton's method !> @param[in] x axial/radial initial guesses for the boundary position !> @param[in] y radial/axial position at which the boundary is searched !> @param[out] xpt position of the boundary found !> @param[inout] direction (1) searches boundary along z, (2) searches boundary along r (0) no boundary found between x(1) and x(2) !--------------------------------------------------------------------------- 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 ! calculates the position of the boundary ( where the weight changes sign ) ! between x(1) and x(2) at 2nd coordinate y xptm=x(1) xpt=x(2) i=0 ! direction=1 finds cross-point along z if(direction .eq.1) Call dom_weight(xptm,y,wold) ! direction=2 finds cross-point along r if(direction .eq.2) Call dom_weight(y,xptm,wold) if(direction .eq.1) Call dom_weight(xpt,y,w) if(direction .eq.2) Call dom_weight(y,xpt,w) ! if the weight doesn't change sign there is no cross point if(w*wold.gt.0) then direction=0 return End If ! Find the cross-point Do while( i .le.100 .and. abs(w).gt.1e-9) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold xptm=xpt-w*(xpt-xptm)/(w-wold) wold=w temp=xptm xptm=xpt xpt=temp i=i+1 if(direction .eq.1) Call dom_weight(xpt,y,w) if(direction .eq.2) Call dom_weight(y,xpt,w) End do if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0 End Subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> calculates the boundary limit between x(1) and x(2) using dichotomy method !> @param[in] x axial/radial initial guesses for the boundary position !> @param[in] y radial/axial position at which the boundary is searched !> @param[out] xpt position of the boundary found !> @param[inout] direction (1) searches boundary along z, (2) searches boundary along r (0) no boundary found between x(1) and x(2) !--------------------------------------------------------------------------- Subroutine Find_crosspointdico(x,y,xpt, direction) ! calculates the boundary limit between x(1) and x(2) using dichotomy method Real(kind=db):: x(2), y Real(kind=db):: xpt Real(kind=db):: w1, w2, w3 Real(kind=db):: x1, x2, x3 Integer, Intent(INOUT):: direction Integer:: i ! calculates the position of the boundary ( where the weight changes sign ) ! between x(1) and x(2) at 2nd coordinate y x1=x(1) x2=x(2) i=0 select case(direction) case(1) ! direction=1 finds cross-point along z Call dom_weight(x1,y,w1) Call dom_weight(x2,y,w2) case(2) ! direction=2 finds cross-point along r Call dom_weight(y,x1,w1) Call dom_weight(y,x2,w2) end select ! if the weight doesn't change sign there is no cross point if(w1*w2.ge.0) then direction=0 xpt=x2 return End If ! Find the cross-point Do while( i .le. 1000 .and. abs((x2-x1)/(x(2)-x(1))).gt.1e-14) x3=0.5*(x1+x2) i=i+1 select case(direction) case(1) ! direction=1 finds cross-point along z Call dom_weight(x3,y,w3) case(2) ! direction=2 finds cross-point along r Call dom_weight(y,x3,w3) end select if(w1*w3.gt.0)then! we are in a region were there is no change of sign x1=x3 w1=w3 else x2=x3 w2=w3 end if !if(abs(w3).lt.1e-14.and. w3 .ge.0)Exit End do xpt=x3 if(xpt .ge. x(2) .or. xpt .le. x(1) ) direction=0 End Subroutine !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> returns the total weight which is the same as the Dirichlet weight for a boundary !> defined with Rvaschev functions !> @param[in] z axial position !> @param[in] r radial position !> @param[out] w dirichlet weight !> @param[out] idwall id of the closest boundary if the particle is outside the domain !--------------------------------------------------------------------------- SUBROUTINE geom_rvaschtot(z,r,w,idwall) ! returns the total weight which is the same as the Dirichlet weight for a boundary ! defined with Rvaschev functions Use splinebound, ONLY: spline_w Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(0:,:) Real(kind=db), allocatable:: w2(:,:) Real(kind=db), allocatable:: w3(:,:) INTEGER, optional, INTENT(OUT):: idwall(:) INTEGER:: sw2 sw2=size(w,2) if(present(idwall)) then idwall=0 allocate(w2(1:size(w,1),1:size(w,2))) allocate(w3(1:size(w,1),1:size(w,2))) call Dirichlet_weight(z,r,w2,w3) ! gives total weight where(w2(1,:).le.0) idwall=1 where(w3(1,:).le.0) idwall=2 call Combine(w2, w3, w, -1) else call Dirichlet_weight(z,r,w) end if End SUBROUTINE geom_rvaschtot Subroutine gstd(z,r,gtilde,w) ! g tilde function added to rhs of poisson solver for the standard coaxial configuration ! for the default weight function Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: belowtmp(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,1)) Real(kind=db):: denom(size(r,1)) if (above1.eq.0) then gtilde=0 gtilde(0,:)=Phiup RETURN end if ! Weight functions necessary for calculation of g ! coaxial insert call cyllweight(z,r,belowtmp, r_a, above1) SELECT CASE (walltype) CASE (0) ! top cylinder call cyllweight(z,r,abovetmp, r_b, above2) CASE DEFAULT ! Ellipse as in gt170 call ellipseweight(z,r,abovetmp,r_0,z_0,invr_z,invr_r,Interior) END SELECT ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(0,:)=(Phidown*abovetmp(0,:) + Phiup*belowtmp(0,:) ) / & & (abovetmp(0,:)+belowtmp(0,:)) If(size(gtilde,1) .gt. 2) then ! first derivative denom=(abovetmp(0,:)+belowtmp(0,:))**2 gtilde(1,:)=(Phiup-Phidown)*(belowtmp(1,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(1,:)) / & & denom ! Axial derivative gtilde(2,:)=(Phiup-Phidown)*(belowtmp(2,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(2,:)) / & & denom ! Radial derivative End If End subroutine SUBROUTINE gUpDown(z,r,gtilde,w) ! g tilde function added to rhs of poisson solver by combining two boundaries set at different potentials Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: belowtmp(0:size(gtilde,1)-1,size(r,1)), abovetmp(0:size(gtilde,1)-1,size(r,1)) Real(kind=db):: denom(size(r,1)) ! Weight functions necessary for calculation of g call Dirichlet_weight(z,r,belowtmp,abovetmp) ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(0,:)=(Phidown*abovetmp(0,:) + Phiup*belowtmp(0,:) ) / & & (abovetmp(0,:)+belowtmp(0,:)) If(size(gtilde,1) .gt. 2) then ! first derivative denom=(abovetmp(0,:)+belowtmp(0,:))**2 gtilde(1,:)=(Phiup-Phidown)*(belowtmp(1,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(1,:)) / & & denom ! Axial derivative gtilde(2,:)=(Phiup-Phidown)*(belowtmp(2,:)*abovetmp(0,:)-belowtmp(0,:)*abovetmp(2,:)) / & & denom ! Radial derivative End If END SUBROUTINE gUpDown Subroutine gtest(z,r,gtilde,w) ! calculates the Poisson gtilde term for testing the solver on a new geometry ! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(0:,:) Real(kind=db):: wtmp(0:size(gtilde,1)-1,size(gtilde,2)) ! !gtilde=0 !return if(present(w)) then wtmp=w else call Dirichlet_weight(z,r,wtmp) end if gtilde(0,:)=(sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))+2) !gtilde(0,:)=1+r**6+2*z**6 !gtilde(0,:)=exp(r*z) If(size(gtilde,1) .gt. 1) then ! first derivative gtilde(1,:)=pi/(test_pars%Lz)*cos(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) gtilde(2,:)=pi/(test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*cos(pi*(r-test_pars%r0)/(test_pars%Lr))*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) ! gtilde(2,:)=6*r**5*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) ! gtilde(1,:)=12*z**5*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) ! !gtilde(2,:)=z*exp(r*z)*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) ! !gtilde(1,:)=r*exp(r*z)*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) End If gtilde(0,:)=gtilde(0,:)*(1-wtmp(0,:)) End subroutine Subroutine manufsolution(z,r,gtilde) ! calculates the Poisson gtilde term for testing the solver on a new geometry ! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(0:,:) ! !gtilde=0 !return gtilde(0,:)=(sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))+2) !gtilde(0,:)=1+r**6+2*z**6 !gtilde(0,:)=exp(r*z) If(size(gtilde,1) .gt. 1) then ! first derivative gtilde(1,:)=-pi/(test_pars%Lz)*cos(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr)) gtilde(2,:)=-pi/(test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*cos(pi*(r-test_pars%r0)/(test_pars%Lr)) ! gtilde(2,:)=6*r**5*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) ! gtilde(1,:)=12*z**5*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) ! !gtilde(2,:)=z*exp(r*z)*(1-wtmp(0,:))-wtmp(2,:)*gtilde(0,:) ! !gtilde(1,:)=r*exp(r*z)*(1-wtmp(0,:))-wtmp(1,:)*gtilde(0,:) End If gtilde(0,:)=gtilde(0,:) End subroutine Subroutine ftest(z,r,f) ! calculates the Poisson source term for testing the solver on a new geometry ! This uses a manufactured solution of the form phi=sin(pi(z-z0)/Lz)sin(pi(r-r0)/Lr)+2 Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT)::f(0:,:) !f(0,:)=1 !return ! f(0,:)=(pi/test_pars%Lz)**2*sin(pi*(z-test_pars%z0)/(test_pars%Lz))*sin(pi*(r-test_pars%r0)/(test_pars%Lr))& & + (pi/test_pars%Lr)*sin(pi*(z-test_pars%z0)/(test_pars%Lz))& & *( -1/r*(cos(pi*(r-test_pars%r0)/(test_pars%Lr))) + (pi/test_pars%Lr)*sin(pi*(r-test_pars%r0)/(test_pars%Lr))) !f(0,:)=-36*r**4-60*z**4 !f(0,:)=-(z*exp(r*z)/r+z**2*exp(r*z)+r**2*exp(r*z)) End Subroutine ftest !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Reduce the Finite element matrix in the LHS from full spline space to reduced web-spline space !> If A is the standard FEM matrix on the full grid space then the web-spline FEM matrix is etilde A etildet !> with etildet the transposed etilde matrix !> @param[in] full full FEM matrix on the weighted b-spline domain !> @param[out] reduced reduced matrix on the web-spline space !--------------------------------------------------------------------------- subroutine Reducematrix(full,reduced) ! use mumps_bsplines type(mumps_mat):: full, reduced, tempmat1, tempmat2 Integer:: fullrank, reducedrank Integer:: i,j, k Real(kind=db):: val logical:: error call to_mat(full) fullrank=full%rank reducedrank=nbreducedspline !WRITE(*,*) "# web-splines", nbreducedspline !WRITE(*,*) "# splines", fullrank call init(reducedrank,2,reduced,nlpos=.false.) call init(fullrank,2,tempmat1) call init(reducedrank,2,tempmat2) tempmat1=mmx_mumps_mat_loc(full,etildet,(/fullrank,fullrank/),(/fullrank,reducedrank/),.false.) tempmat2=mmx_mumps_mat_loc(etildet,tempmat1,(/fullrank,reducedrank/),(/fullrank,reducedrank/),.true.) do i=1,reducedrank do k=tempmat2%irow(i),tempmat2%irow(i+1)-1 j=tempmat2%cols(k) call putele(reduced, i, j, tempmat2%val(k)) end do end do call to_mat(reduced) end subroutine !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !--------------------------------------------------------------------------- !> @author !> Guillaume Le Bars EPFL/SPC ! ! DESCRIPTION: !> !> @brief !> Com putes the sparse matrix matrix multiplication using MKL A B=C !> @param[in] mata sparse matrix A in csr representation !> @param[in] matb sparse matrix B in csr representation !> @param[in] ranka Rank of matrix A !> @param[in] rankb Rank of matrix B !> @param[in] transpose Is A transposed !> @param[out] matc sparse matrix C in csr representation !--------------------------------------------------------------------------- FUNCTION mmx_mumps_mat_loc(mata, matb, ranka, rankb, transpose) RESULT(matc) ! ! Return product matc=mata*matb ! All matrices are represented in csr sparse representation ! Use mumps_bsplines TYPE(mumps_mat) :: mata,matb,matc INTEGER:: ranka(2), rankb(2) ! nb of (rows,columns) for each matrix INTEGER :: n, info, m, request, sort LOGICAL:: transpose Character(1):: T ! m=ranka(1) n=ranka(2) T='N' if (transpose)then T='T' m=ranka(2) n=ranka(1) end if !#ifdef MKL ! if(associated(matc%val)) deallocate(matc%val) if(associated(matc%cols)) deallocate(matc%cols) if(associated(matc%irow)) deallocate(matc%irow) allocate(matc%irow(m+1)) allocate(matc%cols(1)) allocate(matc%val(1)) ! Compute the size of the resulting matrix and the indices of the non-zero elements request=1 sort=7 CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow, & & matb%val, matb%cols, matb%irow, & & matc%val, matc%cols, matc%irow, & & 1, info) if(info .gt. 0) WRITE(*,'(a,i4)') " Error in mmx_mumps_mat_loc: ", info if(associated(matc%val)) deallocate(matc%val) if(associated(matc%cols)) deallocate(matc%cols) allocate(matc%val(matc%irow(m+1)-1)) allocate(matc%cols(matc%irow(m+1)-1)) matc%val=0 request=2 ! Do the actual matrix matrix multiplication CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow, & & matb%val, matb%cols, matb%irow, & & matc%val, matc%cols, matc%irow, & & 1, info) if(info .ne. 0) WRITE(*,'(a,i4)') " Error in mmx_mumps_mat_loc: ", info if (transpose)then matc%rank=ranka(2) else matc%rank=ranka(1) end if ! END FUNCTION mmx_mumps_mat_loc END MODULE geometry