diff --git a/src/geometry_mod.f90 b/src/geometry_mod.f90 index e29c8e3..37faee3 100644 --- a/src/geometry_mod.f90 +++ b/src/geometry_mod.f90 @@ -1,1209 +1,1207 @@ !------------------------------------------------------------------------------ ! 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 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 real(kind=db):: z0 real(kind=db):: r0 real(kind=db):: Lz real(kind=db):: Lr end type Integer, save:: testkr=1, testkz=1 - Logical, save:: nlweb=.false. + Logical, save:: nlweb=.true. 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 Integer, save, Allocatable:: linkedspline(:,:) ! Array containing the lowerleft linked spline in case of boundary spline Real(kind=db), save, allocatable:: gridwgeom(:,:) ! Stores the 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 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(test_params), save:: test_pars PROCEDURE(geom_eval), POINTER:: dirichlet_weight => NULL() PROCEDURE(geomtot_eval), POINTER:: domain_weight => NULL() PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() 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, L_r, L_z, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight, testkr, testkz + 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 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, 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 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 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 (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 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 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 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) end if ALLOCATE(gridwgeom((nrz(1)+1)*(nrz(2)+1),0:2)) gridwgeom=0 ALLOCATE(gridwdom((nrz(1)+1)*(nrz(2)+1),0:2)) gridwdom=0 ALLOCATE(gtilde((nrz(1)+1)*(nrz(2)+1),0:2)) gtilde=0 Call geom_weight (vec1, vec2, gridwgeom) CALL total_gtilde(vec1, vec2, gtilde, gridwgeom) Call dom_weight (vec1, vec2, gridwdom) end Subroutine geom_init Subroutine geom_diag(File_handle, str, rnorm) use mpi Use futils use weighttypes 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_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',gridwdom) Call putarr(File_handle, trim(grpname)//'/dirichletweight',gridwgeom) Call putarr(File_handle, trim(grpname)//'/gtilde',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) END IF End subroutine geom_diag 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. 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 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 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 !!$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 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*(jcellr-icellr-m)/(n-m) 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*(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. 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 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 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(1:1,0:size(w,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 !--------------------------------------------------------------------------- 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(1:1,0:size(w,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(size(w,1),1:1) INTEGER, optional, INTENT(OUT):: idwall(:) call domain_weight(z,r,wtmp,idwall=idwall) 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 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 (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 SUBROUTINE classification(z, r, ctype) ! classify 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 dom_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) ! calculates the gauss integration points for the FEM method for any cell type ! takes care of boundary cells as well 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,wlu,wld, 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 ! 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=4 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),min(xptdown,xptup),max(xptdown,xptup), 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 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 dom_weight(zg(1),rgrid(j),wld) rmin=rgrid(j) if (wld .le. 0) rmin = rgrid(j+1) Do k=1,ngauss(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 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) ! calculates the boundary limit between x(1) and x(2) using Newton's method Real(kind=db):: x(2), y Real(kind=db):: xptm, xpt, temp Real(kind=db):: w, wold 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 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 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.gt.0) then direction=0 xpt=x2 return End If ! Find the cross-point Do while( i .le. 500 .and. abs(x2-x1).gt.1e-13) 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 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(size(w,1),1:size(w,2))) allocate(w3(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 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(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 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(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 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,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,w) ! calculates the Poisson gtilde term for testing the solver on a new geometry 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(size(gtilde,1),0:size(gtilde,2)-1) ! 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) If(size(gtilde,2) .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) End If gtilde(:,0)=gtilde(:,0)*(1-wtmp(:,0)) End subroutine Subroutine ftest(z,r,f) ! calculates the Poisson source term for testing the solver on a new geometry Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT)::f(:,0:) ! 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))) End Subroutine ftest subroutine Reducematrix(full,reduced) ! Reduce the Finite element matrix from full spline space to reduced web-spline space use mumps_bsplines type(mumps_mat):: full, reduced, tempmat1, tempmat2 Integer:: fullrank, reducedrank Integer:: i,j, k call to_mat(full) fullrank=full%rank reducedrank=nbreducedspline call init(reducedrank,2,reduced) 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 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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)) request=1 sort=7 CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow(1), & & matb%val, matb%cols, matb%irow(1), & & matc%val, matc%cols, matc%irow(1), & & 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)) request=2 CALL mkl_dcsrmultcsr(T, request, sort, ranka(1), ranka(2), rankb(2), & & mata%val, mata%cols, mata%irow(1), & & matb%val, matb%cols, matb%irow(1), & & matc%val, matc%cols, matc%irow(1), & & 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