diff --git a/.gitignore b/.gitignore index da52388..2b10b2a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,13 @@ *.h5 *.h5_old *_genmod.f90 *.mod *.o *.out *.m~ -dataanalysis/**/* \ No newline at end of file +*.trace +*.log +dataanalysis/**/* +src/release +src/debug +wk/advi diff --git a/src/geometry_mod.f90 b/src/geometry_mod.f90 index a36c646..92eb9be 100644 --- a/src/geometry_mod.f90 +++ b/src/geometry_mod.f90 @@ -1,840 +1,843 @@ !------------------------------------------------------------------------------ ! 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 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),ALLOCATABLE, SAVE ::IVand (:,:) 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 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 PUBLIC:: geom_weight INTERFACE geom_weight MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2 END INTERFACE geom_weight - NAMELIST /geom_params/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb + NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb contains SUBROUTINE read_geom(Fileid, rnorm, rgrid, Potinn, Potout) + Use mpi + Real(kind=db):: rnorm, rgrid(:), Potinn, Potout - Integer:: Fileid + Integer:: Fileid, mpirank, ierr + CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) Rewind(Fileid) READ(Fileid, geom_params) - WRITE(*, geom_params) + if(mpirank .eq. 0) WRITE(*, geom_params) !! Normalizations and initialization of geometric variables r_a=r_a/rnorm r_b=r_b/rnorm if(r_a .eq. 0 .and. r_b .eq.0) then !! in case no geom_params have been defined r_a=rgrid(1) r_b=rgrid(size(rgrid,1)) 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 Interior=-1 above1=1 above2=-1 Phidown=Potinn Phiup=Potout L_r=L_r/rnorm L_z=L_z/rnorm 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 CompgBoundary(vec1,vec2,gtilde) Call geom_weight(vec1, vec2, gridwgeom) end Subroutine init_geom Subroutine diag_geom(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), "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 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, eij, indexdistance, eij1 integer:: l,m, n, linkedi type(innerspline):: innersplinelist real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:) real(kind=db),allocatable :: etildemat(:,:) !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(nbinspline,nrank(1)*nrank(2))) !etilde=0 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 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 linked splines 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) if(bsplinetype(j) .ne. 0) cycle d=0 where(Ibsplinetype) 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=sqrt(real((jcellz-icellz)**2 + (jcellr-icellr)**2,kind=db)) if(indexdistance .gt. 20)then Write(*,*) 'Error on system conditioning, the number of radial or axial points is too low' stop end if do n=0,norder(2) do i=0,norder(1) eij1=1 !eij=1 do m=0,norder(2) if(n.eq.m) cycle eij1=eij1*(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 eij1=eij1*(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 !findloc(innersplinelist%mu,k,1) eij1=eij1*innersplinelist%weight(k+i+n*nrank(1)) !eij=eij*innersplinelist%weight(k+i+n*nrank(1)) call putele(etilde,linkedi,j,eij1) call putele(etildet,j,linkedi,eij1) 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,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 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), norder(2) Integer:: i, j, mu, imin, imax, jmin, jmax Integer, Allocatable:: splinespan(:,:) call get_dim(spl2, nrank, nrz, norder) bsptype=0 Allocate(splinespan(norder(1)+1,norder(2)+1)) Do mu=1,nrank(1)*nrank(2) i=mod(mu-1,nrank(1))+1 j=(mu-1)/nrank(1)+1 splinespan=-1 imin=max(1,i-norder(1)) imax=min(nrz(1),i) jmin=max(1,j-norder(2)) jmax=min(nrz(2),j) splinespan(1:(imax-imin+1),1:(jmax-jmin+1))=gridctype(imin:imax,jmin:jmax) if ( ANY( splinespan==1 ) ) then bsptype(mu)=1 else if (sum(splinespan).eq. -(norder(1)+1)*(1+norder(2))) then 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), walltmp(1:1,1:1), elliptmp(1:1,1:1) ztmp=z rtmp=r call wallweight(ztmp,rtmp,walltmp, r_a, above1) if (walltype .ne. 0) then call ellipseweight(ztmp,rtmp,elliptmp,Interior) else call wallweight(ztmp,rtmp,elliptmp, r_b, above2) end if if(interior.eq.0 .and. above2.eq.0) then w=walltmp(1,1) else if (above1 .eq. 0) then w=elliptmp(1,1) else w=walltmp(1,1)+elliptmp(1,1)-sqrt(walltmp(1,1)**2+elliptmp(1,1)**2) ! weight at position r,z end if 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), walltmp(1:1,0:size(w,1)-1), elliptmp(1:1,0:size(w,1)-1) Real(kind=db):: squareroot, denom ztmp=z rtmp=r call wallweight(ztmp,rtmp,walltmp, r_a, above1) if (walltype .ne. 0) then call ellipseweight(ztmp,rtmp,elliptmp,Interior) else call wallweight(ztmp,rtmp,elliptmp, r_b, above2) end if if(interior.eq.0 .and. above2.eq.0) then w=walltmp(1,:) else if (above1 .eq. 0) then w=elliptmp(1,:) else denom=walltmp(1,0)**2+elliptmp(1,0)**2 squareroot=sqrt(denom) w(0)=walltmp(1,0)+elliptmp(1,0)-squareroot ! weight at position r,z If(size(w,1) .gt. 1) then ! first derivative w(1)=walltmp(1,1)+elliptmp(1,1)-(elliptmp(1,1)*elliptmp(1,0)+walltmp(1,1)*walltmp(1,0))/& & squareroot ! z derivative of w w(2)=walltmp(1,2)+elliptmp(1,2)-(elliptmp(1,2)*elliptmp(1,0)+walltmp(1,2)*walltmp(1,0))/& & squareroot ! r derivative of w End If If (size(w,1).gt.3) then ! second derivative w(3)=walltmp(1,3)+elliptmp(1,3)-((elliptmp(1,1)*elliptmp(1,0)-walltmp(1,1)*walltmp(1,0))**2+& & (elliptmp(1,3)*elliptmp(1,0)+walltmp(1,3)*walltmp(1,0))*denom)/& & (squareroot*denom)! zz derivative of w w(4)=walltmp(1,4)+elliptmp(1,4)-((elliptmp(1,4)*elliptmp(1,0)-walltmp(1,4)*walltmp(1,0))*denom+& & elliptmp(1,1)*elliptmp(1,2)*walltmp(1,0)**2 + walltmp(1,1)*walltmp(1,2)*elliptmp(1,0)**2- & & elliptmp(1,1)*elliptmp(1,0)*walltmp(1,2)*walltmp(1,0) - & & elliptmp(1,2)*elliptmp(1,0)*walltmp(1,1)*walltmp(1,0) )/& & (squareroot*denom)! zr derivative of w w(5)=walltmp(1,5)+elliptmp(1,5)-((elliptmp(1,2)*elliptmp(1,0)-walltmp(1,2)*walltmp(1,0))**2+& & (elliptmp(1,5)*elliptmp(1,0)+walltmp(1,5)*walltmp(1,0))*denom)/& & (squareroot*denom)! rr derivative of w End If End if End SUBROUTINE geom_weight1 SUBROUTINE geom_weight2(z,r,w) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,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 wallweight(z,r,walltmp, r_a, above1) if (walltype .ne. 0) then call ellipseweight(z,r,elliptmp,Interior) else call wallweight(z,r,elliptmp, r_b, above2) 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 If (size(w,2).gt.3) then ! second derivative w(:,3)=walltmp(:,3)+elliptmp(:,3)-((elliptmp(:,1)*elliptmp(:,0)-walltmp(:,1)*walltmp(:,0))**2+& & (elliptmp(:,3)*elliptmp(:,0)+walltmp(:,3)*walltmp(:,0))*denom)/& & (squareroot*denom)! zz derivative of w w(:,4)=walltmp(:,4)+elliptmp(:,4)-((elliptmp(:,4)*elliptmp(:,0)-walltmp(:,4)*walltmp(:,0))*denom+& & elliptmp(:,1)*elliptmp(:,2)*walltmp(:,0)**2 + walltmp(:,1)*walltmp(:,2)*elliptmp(:,0)**2- & & elliptmp(:,1)*elliptmp(:,0)*walltmp(:,2)*walltmp(:,0) - & & elliptmp(:,2)*elliptmp(:,0)*walltmp(:,1)*walltmp(:,0) )/& & (squareroot*denom)! zr derivative of w w(:,5)=walltmp(:,5)+elliptmp(:,5)-((elliptmp(:,2)*elliptmp(:,0)-walltmp(:,2)*walltmp(:,0))**2+& & (elliptmp(:,5)*elliptmp(:,0)+walltmp(:,5)*walltmp(:,0))*denom)/& & (squareroot*denom)! rr derivative of w End If End if End SUBROUTINE geom_weight2 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) 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) .lt. 0 ) then ctype=0 return End If If(weights(1,1)*weights(3,1) .lt. 0 ) then ctype=0 return End If If(weights(2,1)*weights(4,1) .lt. 0 ) then ctype=0 return End If If(weights(3,1)*weights(4,1) .lt. 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=.false. 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 hasmaxpoint=.False. 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 wallweight(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 If (size(w,2).gt.3) then ! second derivative w(:,3)= 0 ! zz derivative w(:,4)= 0 ! zr derivative w(:,5)= 0 ! rr derivative End If End subroutine wallweight Subroutine ellipseweight(z,r,w, Interior) Real(kind=db):: z(:),r(:),w(:,0:) Integer:: Interior w(:,0)=Interior*(1-(r-r_0)**2*invr_r**2-(z-z_0)**2*invr_z**2) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=-2*(z-z_0)*invr_z**2*Interior ! z derivative of w w(:,2)=-2*(r-r_0)*invr_r**2*Interior ! r derivative of w End If If (size(w,2).gt.3) then ! second derivative w(:,3)=-2*invr_z**2*Interior ! zz derivative w(:,4)=0 ! zr derivative w(:,5)=-2*invr_r**2*Interior ! rr derivative End If End subroutine ellipseweight Subroutine CompgBoundary(z,r,gtilde) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) if (walltype .eq. -1) then call gtest(z,r,gtilde) else call gstd(z,r,gtilde) end if End subroutine Subroutine gstd(z,r,gtilde) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: walltmp(size(r,1),0:size(gtilde,2)-1), elliptmp(size(r,1),0:size(gtilde,2)-1) Real(kind=db):: denom(size(r,1)) call wallweight(z,r,walltmp, r_a, above1) if (walltype .eq. 1) then call ellipseweight(z,r,elliptmp,Interior) else call wallweight(z,r,elliptmp, r_b, above2) end if ! Extension to the whole domain of the boundary conditions gtilde(:,0)=(Phidown*elliptmp(:,0) + Phiup*walltmp(:,0) ) / & & (elliptmp(:,0)+walltmp(:,0)) If(size(gtilde,2) .gt. 1) then ! first derivative denom=(elliptmp(:,0)+walltmp(:,0))**2 gtilde(:,1)=(Phiup-Phidown)*(walltmp(:,1)*elliptmp(:,0)-walltmp(:,0)*elliptmp(:,1)) / & & denom ! Axial derivative gtilde(:,2)=(Phiup-Phidown)*(walltmp(:,2)*elliptmp(:,0)-walltmp(:,0)*elliptmp(:,2)) / & & denom ! Radial derivative End If End subroutine Subroutine gtest(z,r,gtilde) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) ! 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) gtilde(:,2)=-2*pi/L_r*sin(2*pi*z/L_z)*sin(2*pi*r/L_r) End If 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(:) Integer:: fullrank, reducedrank Integer:: i fullrank=full%rank reducedrank=nbreducedspline call init(reducedrank,2,reduced) call init(fullrank,2,tempmat) allocate(temparray(fullrank)) !$OMP parallel private(i), private(temparray) !$OMP DO do i=1,fullrank call getcol(full,i,temparray) temparray=vmx(etilde,temparray) !$OMP critical (tempmat_updt) call putcol(tempmat,i,temparray) !$OMP end critical (tempmat_updt) end do !$OMP end do !$OMP end parallel call to_mat(tempmat) !$OMP parallel private(i), private(temparray) !allocate(temparray(nrank(1)*nrank(2)),tempresult(nrank(1)*nrank(2))) !$OMP DO do i=1,reducedrank call getcol(etildet,i,temparray) temparray=vmx(tempmat,temparray) !$OMP critical (reduccedmat_updt) call putcol(reduced,i,temparray(1:reducedrank)) !$OMP end critical (reduccedmat_updt) end do !$OMP end do !$OMP end parallel end subroutine END MODULE geometry \ No newline at end of file