diff --git a/src/geometry_mod.f90 b/src/geometry_mod.f90 index 134326c..3d95594 100644 --- a/src/geometry_mod.f90 +++ b/src/geometry_mod.f90 @@ -1,1282 +1,1338 @@ !------------------------------------------------------------------------------ ! 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),save:: z_0=0, r_0=0, invr_r=0, invr_z=0, r_a=0, r_b=0, z_r=1, r_r=1, alpha=0 Real(kind=db),save:: r_bLeft=0, r_bRight=0 Real(kind=db), save:: Phidown=0,Phiup=0 Real(kind=db), save:: L_r=0.16, L_z=0.24 Logical, save:: nlweb=.false. Integer, save::Interior=-1 Integer, save:: above1=1, above2=-1 Integer, save :: walltype = 0 Integer, save, Allocatable:: bsplinetype(:) ! Array containing the inner/outer type for each bspline Integer, save, Allocatable:: gridcelltype(:,:) ! Array containing the inner/outer type for each gridcell Real(kind=db), save, allocatable:: gridwgeom(:,:) ! Stores the geometric weight at the grid points Real(kind=db), save, allocatable:: gtilde(:,:) ! Stores the extension to the domain of the boundary conditions type(mumps_mat):: etilde ! Matrix of extendend web splines type(mumps_mat):: etildet ! Transpose of Matrix of extendend web splines integer,save :: nbreducedspline ! Number of splines in the reduced set PROCEDURE(geom_eval), POINTER:: total_weight => NULL() PROCEDURE(gtilde_eval), POINTER:: total_gtilde => NULL() PUBLIC:: geom_weight ABSTRACT INTERFACE SUBROUTINE gtilde_eval(z,r,g) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(:,0:) END SUBROUTINE SUBROUTINE geom_eval(z,r,w,wupper) USE constants Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), INTENT(OUT), OPTIONAL:: wupper(:,0:) END SUBROUTINE END INTERFACE INTERFACE geom_weight MODULE PROCEDURE geom_weight0, geom_weight1, geom_weight2 END INTERFACE geom_weight NAMELIST /geomparams/ z_0, r_0, z_r, r_r, r_a, r_b, walltype, L_r, L_z, nlweb, Interior, above1, above2, alpha, r_bLeft, r_bRight contains SUBROUTINE read_geom(Fileid, rnorm, rgrid, Potinn, Potout) Use mpi Real(kind=db):: rnorm, rgrid(:), Potinn, Potout Integer:: Fileid, mpirank, ierr, istat character(len=1000) :: line !type(spline_boundary):: border0 !Real(kind=db):: cpoints(2,3),t(1), distance, pos !Real(kind=db), ALLOCATABLE::result(:,:,:) CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) Rewind(Fileid) READ(Fileid, geomparams, iostat=istat) if (istat.gt.0) then if(mpirank .eq. 0) then backspace(Fileid) read(Fileid,fmt='(A)') line write(*,'(A)') & 'Invalid line in geomparams: '//trim(line) end if call MPI_Abort(MPI_COMM_WORLD, -1, ierr) stop end if if(mpirank .eq. 0) WRITE(*, geomparams) !! Normalizations and initialization of geometric variables r_a=r_a/rnorm r_b=r_b/rnorm r_bLeft=r_bLeft/rnorm r_bRight=r_bRight/rnorm if(r_a .eq. 0 .and. r_b .eq.0) then !! in case no geom_params have been defined r_a=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 Phidown=Potinn Phiup=Potout L_r=L_r/rnorm L_z=L_z/rnorm SELECT CASE (walltype) CASE(-1) ! test case total_gtilde=>gtest total_weight=>geom_weightstd CASE (2) ! coaxial cylinder and top ellipse with cylinder extensions total_gtilde=>gUpDown total_weight=>geom_w2 CASE (3) ! Two ellipses with "parallel" tangents with cylinders total_gtilde=>gUpDown total_weight=>geom_w3 CASE (4) ! Two ellipses with same radii with cylinders total_gtilde=>gUpDown total_weight=>geom_w4 CASE (5) ! Two ellipses with same radii with cylinders and total_gtilde=>gUpDown total_weight=>geom_w5 CASE (6) - ! Two ellipses with same radii with cylinders and + ! circular coaxial tilted ellipse right Dirichlet total_gtilde=>gUpDown total_weight=>geom_w6 + CASE (7) + ! circular coaxial tilted ellipse right and left Dirichlet + total_gtilde=>gUpDown + total_weight=>geom_w7 CASE DEFAULT ! Ellipse as in gt170 standard weight and straight coaxial configs total_gtilde=>gstd total_weight=>geom_weightstd END SELECT !cpoints=reshape(real([0.0, 1.0, 0.0, 2.0, 0.0, 3.0],kind=db),(/2,3/)) !call setspline_boundary(border0, cpoints, 1, 0.0_db) !t(1)=0.75 !call eval(border0, t, result) !call dist(border0,real((/1.0,1.5/),kind=db),distance,pos) end subroutine read_geom SUBROUTINE init_geom(spl2, vec1, vec2) type(spline2d):: spl2 real(kind=db):: vec1(:),vec2(:) Real(kind=db), Allocatable:: zgrid(:),rgrid(:) Integer:: nrank(2), nrz(2) Call get_dim(spl2, nrank, nrz) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) ! create a table of the calssification for each cell Allocate(gridcelltype(nrz(1),nrz(2))) Call classifycell(zgrid, rgrid, gridcelltype) if (nlweb) then Allocate(bsplinetype(nrank(1)*nrank(2))) Call classifyspline(spl2, gridcelltype, bsplinetype) Call buildetilde(spl2,bsplinetype,gridcelltype,etilde) end if ALLOCATE(gridwgeom((nrz(1)+1)*(nrz(2)+1),0:2)) gridwgeom=0 ALLOCATE(gtilde((nrz(1)+1)*(nrz(2)+1),0:2)) gtilde=0 CALL total_gtilde(vec1, vec2, gtilde) Call geom_weight (vec1, vec2, gridwgeom) end Subroutine init_geom Subroutine geom_diag(File_handle, str, rnorm) Use mpi Use futils Integer:: File_handle Real(kind=db):: rnorm Character(len=*):: str CHARACTER(len=128):: grpname Integer:: ierr, mpirank CALL MPI_COMM_RANK(MPI_COMM_WORLD, mpirank, ierr) IF(mpirank .eq. 0) THEN Write(grpname,'(a,a)') trim(str),"/geometry" If(.not. isgroup(File_handle, trim(grpname))) THEN CALL creatg(File_handle, trim(grpname)) END IF Call attach(File_handle, trim(grpname), "r_a", r_a*RNORM) Call attach(File_handle, trim(grpname), "r_b", r_b*RNORM) Call attach(File_handle, trim(grpname), "z_0", z_0*RNORM) Call attach(File_handle, trim(grpname), "r_0", r_0*RNORM) Call attach(File_handle, trim(grpname), "r_r", r_r*RNORM) Call attach(File_handle, trim(grpname), "z_r", z_r*RNORM) Call attach(File_handle, trim(grpname), "interior", interior) Call attach(File_handle, trim(grpname), "above1", above1) Call attach(File_handle, trim(grpname), "above2", above2) Call attach(File_handle, trim(grpname), "walltype", walltype) Call putarr(File_handle, trim(grpname)//'/geomweight',gridwgeom) END IF End subroutine geom_diag Subroutine buildetilde(spl2, bsplinetype, celltype, etilde) USE mumps_bsplines type(spline2d):: spl2 Integer:: bsplinetype(:) Integer:: celltype(:,:) type(mumps_mat):: etilde Logical, Allocatable:: Ibsplinetype(:) Integer:: i,j,k, icellz, icellr, jcellz, jcellr, nrank(2), nrz(2), norder(2) real(kind=db), allocatable:: zgrid(:), rgrid(:) real(kind=db):: wgeomi, indexdistance, eij integer:: l,m, n, linkedi type(innerspline):: innersplinelist real(kind=db), allocatable:: rgridmesh(:),zgridmesh(:), d(:), hz(:), cz(:), hr(:), cr(:), hmesh(:) !if(allocated(etilde)) deallocate(etilde) nbreducedspline=count(bsplinetype .eq. 1) Call get_dim(spl2,nrank,nrz,norder) ! Obtain grid data drom the spline structure Allocate(zgrid(1:nrz(1)+1)) Allocate(rgrid(1:nrz(2)+1)) zgrid=spl2%sp1%knots(0:nrz(1)) rgrid=spl2%sp2%knots(0:nrz(2)) allocate(innersplinelist%k(nrank(1)*nrank(2)),innersplinelist%weight(nrank(1)*nrank(2))) allocate(Ibsplinetype(nrank(1)*nrank(2))) allocate(rgridmesh(nrank(1)*nrank(2)), zgridmesh(nrank(1)*nrank(2))) allocate(hmesh(nrank(1)*nrank(2))) allocate(d(size(bsplinetype))) Ibsplinetype=.False. call calcsplinecenters(spl2%sp1,cz,hz) call calcsplinecenters(spl2%sp2,cr,hr) Do i=0,nrank(2)-1 zgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cz rgridmesh(i*nrank(1)+1:(i+1)*nrank(1))=cr(i+1) hmesh(i*nrank(1)+1:(i+1)*nrank(1))=hr(i+1)*hz End do ! allocate etilde and its transpose call init(nrank(1)*nrank(2),nbreducedspline,etilde) call init(nrank(1)*nrank(2),nbreducedspline,etildet) k=1 Do i=1,nrank(1)*nrank(2) if(bsplinetype(i) .ne. 1) cycle ! span of this spline is completely outside D ! one cell of bspline i is completely in D icellz=mod(i-1,nrank(1))+1 icellr=(i-1)/(nrank(1))+1 outer: do l=max(1,icellz-norder(1)),min(nrz(1),icellz) do m=max(1,icellr-norder(2)),min(nrz(2),icellr) if (celltype(l,m).eq.1) then call geom_weight((zgrid(l)+zgrid(l+1))/2,(rgrid(m)+rgrid(m+1))/2,wgeomi) EXIT outer end if end do end do outer call putele(etilde,k,i,1/wgeomi) call putele(etildet,i,k,1/wgeomi) innersplinelist%k(i)=k innersplinelist%weight(i)=1/wgeomi k=k+1 !Ibsplinetype(i)=icellz-norder(1) .gt. 0 .and. icellr-norder(2) .gt. 0 .and. icellz .le. nrz(1) .and. icellr .le. nrz(2) Ibsplinetype(i)=icellz+norder(1) .le. nrank(1) .and. icellr+norder(2) .le. nrank(2) if(.not. Ibsplinetype(i)) cycle do m=0,norder(2) do l=0,norder(1) ! Check if all 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) ! find the closest b-spline fully in D for the web method 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. (1+norder(1))*(norder(2)+1))then Write(*,*) 'Error on system conditioning, the number of radial or axial points is too low!' Write(*,*) 'Distance found: ', indexdistance !stop end if do n=0,norder(2) do i=0,norder(1) eij=1 !eij=1 do m=0,norder(2) if(n.eq.m) cycle eij=eij*(rgridmesh(j)-rgridmesh(k+m*nrank(1)))/(rgridmesh(k+n*nrank(1))-rgridmesh(k+m*nrank(1))) end do do l=0,norder(1) if( i .eq. l ) cycle eij=eij*(zgridmesh(j)-zgridmesh(k+l))/(zgridmesh(k+i)-zgridmesh(k+l)) !eij=eij*(jcellr-icellr-m)/(n-m)*(jcellz-icellz-l)/(i-l) end do linkedi=innersplinelist%k(k+i+n*nrank(1)) ! equivalent to findloc, necessary for ifort 17 !findloc(innersplinelist%mu,k,1) eij=eij*innersplinelist%weight(k+i+n*nrank(1)) call putele(etilde,linkedi,j,eij) call putele(etildet,j,linkedi,eij) end do end do end do !!$OMP End parallel do call to_mat(etilde) call to_mat(etildet) end subroutine Subroutine calcsplinecenters(spl,ctrs,heights) use bsplines type(spline1d):: spl real(kind=db), allocatable:: ctrs(:), heights(:) integer:: nrank, nx, order, i, left1, j, left2, left3 real(kind=db):: x1, x2, x3 real(kind=db), allocatable:: fun1(:,:), fun2(:,:), fun3(:,:) real(kind=db),allocatable:: init_heights(:) call get_dim(spl, nrank, nx, order) if (allocated(ctrs)) deallocate(ctrs) if (allocated(heights)) deallocate(heights) allocate(heights(nrank)) allocate(ctrs(nrank)) allocate(fun1(1:order+1,0:1), fun2(1:order+1,0:1), fun3(1:order+1,0:1)) allocate(init_heights(nrank)) init_heights=1.0 ctrs(:)=spl%knots(0:nrank-1) call gridval(spl, ctrs, heights, 0, init_heights) if (order .gt.1) then do i=2,nrank-1 x1=spl%knots(i-order)+1e5*Epsilon(x1) x2=spl%knots(i-1)-1e5*Epsilon(x2) left1=max(0,i-order) left2=min(nx-2,i-2) call basfun(x1,spl,fun1, left1+1) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold call basfun(x2,spl,fun2, left2+1) fun3=1 j=0 Do while( j .le.100) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold x3=(x1+x2)/2 call locintv(spl, x3, left3) call basfun(x3,spl,fun3, left3+1) if( (x2-x1)/2.lt.1e-12) exit if(abs(fun3(i-left3,1)).lt.1e-15) exit if(fun3(i-left3,1)*fun1(i-left1,1).le.0) then fun2=fun3 x2=x3 left2=left3 else fun1=fun3 x1=x3 left1=left3 end if j=j+1 End do ctrs(i)=x3 heights(i)=fun3(i-left3,0) end do end if end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Subroutine classifycell(zgrid,rgrid,gridctype) Real(kind=db):: zgrid(:),rgrid(:) Integer:: gridctype(:,:) Integer::i,j gridctype=-1 !$OMP parallel do private(i) Do j=1,size(rgrid,1)-1 Do i=1,size(zgrid,1)-1 ! Determines the type inner/boundary/outer for each cell Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),& & gridctype(i,j)) End Do End do !$OMP end parallel do End subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine classifyspline(spl2,gridctype,bsptype) type(spline2d):: spl2 Integer:: gridctype(:,:) Integer:: bsptype(:) Integer:: nrank(2), nrz(2), norder(2) Integer:: i, j, mu, imin, imax, jmin, jmax Integer, Allocatable:: splinespan(:,:) call get_dim(spl2, nrank, nrz, norder) ! by default, all splines have part of their support outside the domain D bsptype=0 Allocate(splinespan(norder(1)+1,norder(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 splinespan=-1 ! find the axial span of this spline in cell indices imin=max(1,i-norder(1)) imax=min(nrz(1),i) ! find the radial span of this spline in cell indices jmin=max(1,j-norder(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 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), wtmp(1:1,1:1) ztmp=z rtmp=r call total_weight(ztmp,rtmp,wtmp) w=wtmp(1,1) End SUBROUTINE geom_weight0 SUBROUTINE geom_weight1(z,r,w) Real(kind=db), INTENT(IN):: r,z Real(kind=db), INTENT(OUT):: w(0:) Real(kind=db):: rtmp(1:1), ztmp(1:1) Real(kind=db):: wtmp(1:1,0:size(w,1)-1) ztmp=z rtmp=r call total_weight(ztmp,rtmp,wtmp) w=wtmp(1,:) End SUBROUTINE geom_weight1 SUBROUTINE geom_weight2(z,r,w) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) call total_weight(z,r,w) ! End SUBROUTINE geom_weight2 SUBROUTINE geom_weightstd(z,r,w,wupper) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) Real(kind=db):: squareroot(size(w,1)), denom(size(w,1)) call cyllweight(z,r,walltmp, r_a, above1) SELECT CASE (walltype) CASE (0) call cyllweight(z,r,elliptmp, r_b, above2) CASE DEFAULT call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior) END SELECT if(interior.eq.0 .and. above2.eq.0) then w=walltmp(:,:) else if (above1 .eq. 0) then w=elliptmp(:,:) else denom=walltmp(:,0)**2+elliptmp(:,0)**2 squareroot=sqrt(denom) w(:,0)=walltmp(:,0)+elliptmp(:,0)-squareroot ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=walltmp(:,1)+elliptmp(:,1)-(elliptmp(:,1)*elliptmp(:,0)+walltmp(:,1)*walltmp(:,0))/& & squareroot ! z derivative of w w(:,2)=walltmp(:,2)+elliptmp(:,2)-(elliptmp(:,2)*elliptmp(:,0)+walltmp(:,2)*walltmp(:,0))/& & squareroot ! r derivative of w End If End if End SUBROUTINE geom_weightstd SUBROUTINE geom_w2(z,r,w,wupper) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) call cyllweight(z,r,walltmp, r_a, above1) call ellipsewall(z,r,elliptmp,r_b,above2,r_0, z_0, invr_z, invr_r, Interior) if( present(wupper) ) then wupper=elliptmp w=walltmp else call Combine(elliptmp, walltmp, w, -1) end if End SUBROUTINE geom_w2 SUBROUTINE geom_w3(z,r,w,wupper) ! Defines the geometric weight for two facing ellipses of parallel boundaries with line prolongations !-----\__/------ ! !---\______/---- Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: belowtmp(size(w,1),0:size(w,2)-1), abovetmp(size(w,1),0:size(w,2)-1) ! Weight functions necessary for calculation of top and bottom boundaries ! Ellipse and bottom cylinder combination call ellipsewall(z,r,belowtmp,r_a,1,r_b, z_0, 1/(z_r+r_b-r_a), 1/(z_r+r_b-r_a), 1) ! Ellipse and top cylinder combination call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1) if( present(wupper) ) then wupper=abovetmp w=belowtmp else call Combine(belowtmp, abovetmp, w, -1) end if End SUBROUTINE geom_w3 SUBROUTINE geom_w4(z,r,w,wupper) ! Defines the geometric weight for two facing ellipses of same dimensions with line prolongations !-----\__/------ ! !-----\__/------ Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: belowtmp(size(w,1),0:size(w,2)-1), abovetmp(size(w,1),0:size(w,2)-1) ! Weight functions necessary for calculation of g ! Ellipse and bottom cylinder combination call ellipsewall(z,r,belowtmp,r_a,1,r_a, z_0, invr_z, invr_r, 1) ! Ellipse and top cylinder combination call ellipsewall(z,r,abovetmp,r_b,-1,r_b, z_0, invr_z, invr_r, -1) if( present(wupper) ) then wupper=abovetmp w=belowtmp else call Combine(belowtmp, abovetmp, w, -1) end if End SUBROUTINE geom_w4 SUBROUTINE geom_w5(z,r,w, wupper) ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w5 SUBROUTINE geom_w6(z,r,w, wupper) use basic, ONLY: zgrid, nz ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,0:) Real(kind=db), OPTIONAL:: wupper(:,0:) Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! tilted ellipse call discweight(z,r, w1, zgrid(nz),-1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w2 w=w1 else ! gives total weight call Combine(w1, w2, w, -1) end if End SUBROUTINE geom_w6 + + SUBROUTINE geom_w7(z,r,w, wupper) + use basic, ONLY: zgrid, nz + ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region + ! and left and right dirichlet boundaries + Real(kind=db), INTENT(IN):: r(:),z(:) + Real(kind=db), INTENT(OUT):: w(:,0:) + Real(kind=db), OPTIONAL:: wupper(:,0:) + Real(kind=db):: w1(size(w,1),0:size(w,2)-1), w2(size(w,1),0:size(w,2)-1), w3(size(w,1),0:size(w,2)-1) + + + ! tilted wall + call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) + ! calculate upper right boundary + call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) + ! Intersection between slanted wall and upper right wall + call Combine(w1,w2, w3, -1) + + ! calculate upper left boundary + call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) + ! Union between previous weight w3 and upper left wall + call Combine(w1,w3, w2, 1) + + ! tilted ellipse + call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) + ! Intersection of ellipse and previous weight + call Combine(w1, w2, w3, Interior) + + ! Right Dirichlet + call discweight(z,r, w1, zgrid(nz),-1) + ! Intersection of ellipse and previous weight + call Combine(w1, w3, w2, -1) + + ! Left Dirichlet + call discweight(z,r, w1, zgrid(0),-1) + ! Intersection of ellipse and previous weight + call Combine(w1, w2, w3, -1) + + ! calculate coaxial insert + call cyllweight(z,r,w1,r_a,above1) + + if( present(wupper) ) then + wupper=w2 + w=w1 + else + ! gives total weight + call Combine(w1, w3, w, -1) + end if + + End SUBROUTINE geom_w7 + 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 Real(kind=db):: xptup, xptdown, w, rmin, rmax type(spline1d):: splz, splr splz=spl2%sp1 splr=spl2%sp2 Allocate(zgrid(1:splz%nints+1)) Allocate(rgrid(1:splr%nints+1)) zgrid=splz%knots(0:splz%nints) rgrid=splr%knots(0:splr%nints) hasmaxpoint=.false. If(allocated(xgauss)) deallocate(xgauss) if(allocated(wgauss)) deallocate(wgauss) !Call classification((/zgrid(i),zgrid(i+1)/),(/rgrid(j),rgrid(j+1)/),ctype) ctype=gridcelltype(i,j) If (ctype .ge. 1) then ! we have a normal internal cell Allocate(xgauss(ngauss(1)*ngauss(2),2)) Allocate(wgauss(ngauss(1)*ngauss(2))) gausssize=ngauss(1)*ngauss(2) !Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration DO k=1,ngauss(2) xgauss((k-1)*ngauss(1)+1:k*ngauss(1),1)=zg xgauss((k-1)*ngauss(1)+1:k*ngauss(1),2)=rg(k) wgauss((k-1)*ngauss(1)+1:k*ngauss(1))=wrg(k)*wzg END DO Else If(ctype.eq.0) then ! we have a boundary cell !Call Boundary_Limits(pts) !Do k=1,size(pts,1) ! hasmaxpoint=hasmaxpoint .or. (pts(k,1) .ge. zgrid(i) .and. pts(k,1) .le. zgrid(i+1) & ! & .and. pts(k,2) .ge. rgrid(j) .and. pts(k,2) .le. rgrid(j+1)) !End do If(.not. hasmaxpoint) Then directiondown=1 directionup=1 Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j),xptdown,directiondown) Call Find_crosspoint((/zgrid(i),zgrid(i+1)/),rgrid(j+1),xptup,directionup) select case ( directionup+directiondown) Case (0) ! The intersections are only on the left and right cell boundaries nbzpoints=2 Allocate(zpoints(nbzpoints)) zpoints=(/zgrid(i),zgrid(i+1)/) Case(1) nbzpoints=3 Allocate(zpoints(nbzpoints)) if(directiondown.eq.1) zpoints=(/zgrid(i), xptdown,zgrid(i+1)/) if(directionup .eq. 1) zpoints=(/zgrid(i), xptup,zgrid(i+1)/) Case(2) nbzpoints=3 Allocate(zpoints(nbzpoints)) call geom_weight(zgrid(i),rgrid(j),w) If(w.ge.0) zpoints=(/zgrid(1),min(xptdown,xptup),max(xptdown,xptup) /) If(w.lt.0) zpoints=(/min(xptdown,xptup),max(xptdown,xptup), zgrid(i+1) /) End select Allocate(xgauss(ngauss(1)*ngauss(2)*(nbzpoints-1),2)) Allocate(wgauss(ngauss(1)*ngauss(2)*(nbzpoints-1))) gausssize=ngauss(1)*ngauss(2)*(nbzpoints-1) ! Compute gauss points Do l=1,nbzpoints-1 !CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) Call gauleg(zpoints(l),zpoints(l+1),zg,wzg,ngauss(1)) ! We test if the lower or upper side is in the domain call geom_weight(zg(1),rgrid(j),w) rmin=rgrid(j) if (w .lt. 0) rmin = rgrid(j+1) Do k=1,ngauss(1) direction=2 Call Find_crosspoint((/rgrid(j),rgrid(j+1)/),zg(k),rmax,direction) ! We compute the radial limits at each z position Call gauleg(min(rmin,rmax),max(rmin,rmax),rg,wrg,ngauss(2)) ! We obtain the gauss w and pos for these boundaries xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),1) = zg(k) xgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1),2) = rg wgauss(k+(l-1)*ngauss(1)*ngauss(2) : l*ngauss(2)*ngauss(1) : ngauss(1)) = wrg*wzg(k) End do End Do End If Else gausssize=1 Allocate(xgauss(1:1,2)) Allocate(wgauss(1:1)) !Computation of gauss weight and position in r and z direction for gaussian integration CALL get_gauss(spl2%sp1, ngauss(1), i, zg, wzg) CALL get_gauss(spl2%sp2, ngauss(2), j, rg, wrg) ! Construction of matrix xgauss and wgauss storing the weight and position for 2d gaussian integration xgauss(1,1)=(zg(1)+zg(ngauss(1)))*0.5 xgauss(1,2)=(rg(1)+rg(ngauss(2)))*0.5 wgauss(1)=0 End If If(PRESENT(celltype)) celltype=ctype End Subroutine calc_gauss Subroutine Find_crosspoint(x,y,xpt, direction) Real(kind=db):: x(2), y Real(kind=db):: xptm, xpt, temp Real(kind=db):: w, wold Integer, Intent(INOUT):: direction Integer:: i xptm=x(1) xpt=x(2) i=0 if(direction .eq.1) Call geom_weight(xptm,y,wold) if(direction .eq.2) Call geom_weight(y,xptm,wold) if(direction .eq.1) Call geom_weight(xpt,y,w) if(direction .eq.2) Call geom_weight(y,xpt,w) if(w*wold.gt.0) then direction=0 return End If Do while( i .le.100 .and. w*wold .lt. 0) !Write(*,*) 'i,xpt,xptm,w,wold,',i,xpt,xptm,w,wold if(direction .eq.1) Call geom_weight(xpt,y,w) if(direction .eq.2) Call geom_weight(y,xpt,w) xptm=xpt-w*(xpt-xptm)/(w-wold) wold=w temp=xptm xptm=xpt xpt=temp i=i+1 End do if(xpt .ge. x(2) .or. xpt .le. x(1)) direction=0 End Subroutine SUBROUTINE cyllweight(z,r,w, r_lim, above) Real(kind=db):: z(:),r(:),w(:,0:), r_lim Integer:: above w(:,0)=above*(r-r_lim) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=0 ! z derivative of w w(:,2)=above ! r derivative of w End If End subroutine cyllweight SUBROUTINE discweight(z,r,w, z_lim, right) Real(kind=db):: z(:),r(:),w(:,0:), z_lim Integer:: right w(:,0)=right*(z-z_lim) ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=right ! z derivative of w w(:,2)=0 ! r derivative of w End If End subroutine discweight SUBROUTINE tiltedwall(z,r,w,r0,z0,alpha,left) Real(kind=db):: z(:),r(:),w(:,0:) Integer:: left Real(kind=db):: slope, r0, z0, alpha slope=tan(alpha) w(:,0)=(r-r0-slope*(z-z0)) If(size(w,2) .gt. 1) then ! first derivative w(:,1)= -slope ! z derivative of w w(:,2)= 1 ! r derivative of w End If w=left*w end SUBROUTINE SUBROUTINE ellipsewall(z,r,w,r_lim,above,r_0, z_0, invr_z, invr_r, Interior) Real(kind=db):: z(:),r(:),w(:,0:) Integer:: Interior, above Real(kind=db):: r_lim, r_0,z_0,invr_z,invr_r Real(kind=db):: walltmp(size(w,1),0:size(w,2)-1), elliptmp(size(w,1),0:size(w,2)-1) call cyllweight(z,r,walltmp, r_lim, above) call ellipseweight(z,r,elliptmp,r_0, z_0, invr_z, invr_r, Interior) call combine( walltmp, elliptmp, w, Interior) END SUBROUTINE ellipsewall Subroutine ellipseweight(z,r,w, r_0, z_0, invr_z, invr_r, Interior) Real(kind=db):: z(:),r(:),w(:,0:) Real(kind=db):: r_0,z_0,invr_z,invr_r Integer:: Interior 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 End subroutine ellipseweight Subroutine tiltedellipse(z,r,w, r_0, z_0, invr_z, invr_r, Interior, alpha) Real(kind=db):: z(:),r(:),w(:,0:) Real(kind=db):: r_0,z_0,invr_z,invr_r, alpha, cosa, sina Real(kind=db):: deltar(size(r,1)), deltaz(size(z,1)) Integer:: Interior cosa=cos(alpha) sina=sin(alpha) deltar=(r-r_0) deltaz=(z-z_0) w(:,0)=1-((deltaz*cosa+deltar*sina)*invr_z)**2-((deltaz*sina-deltar*cosa)*invr_r)**2 ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=-2*cosa*(deltaz*cosa+deltar*sina)*invr_z**2-2*sina*(deltaz*sina-deltar*cosa)*invr_r**2 ! z derivative of w w(:,2)=-2*sina*(deltaz*cosa+deltar*sina)*invr_z**2+2*cosa*(deltaz*sina-deltar*cosa)*invr_r**2 ! r derivative of w End If w=w*Interior End subroutine tiltedellipse SUBROUTINE combine(w1, w2, w, union) ! Combines two weights w1 and w2 into w ! If union is 1 this defines the union of the geometric domains ! If union is -1 this defines the intersection of the geometric domains Real(kind=db):: w1(:,0:), w2(:,0:), w(:,0:) Integer:: union ! if 1 defines the union ow weights, if -1 defines intersection Real(kind=db):: squareroot(size(w,1)), denom(size(w,1)) denom=w1(:,0)**2+w2(:,0)**2 squareroot=union*sqrt(denom) w(:,0)=w1(:,0)+w2(:,0)+squareroot ! weight at position r,z If(size(w,2) .gt. 1) then ! first derivative w(:,1)=w1(:,1)+w2(:,1)+(w1(:,1)*w1(:,0)+w2(:,1)*w2(:,0))/squareroot ! z derivative of w w(:,2)=w1(:,2)+w2(:,2)+(w1(:,2)*w1(:,0)+w2(:,2)*w2(:,0))/squareroot ! r derivative of w End If END SUBROUTINE Subroutine gstd(z,r,gtilde) ! standard g function added to rhs of poisson solver Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1) Real(kind=db):: denom(size(r,1)) if (above1.eq.0) then gtilde=0 gtilde(:,0)=Phiup RETURN end if ! Weight functions necessary for calculation of g ! coaxial insert call cyllweight(z,r,belowtmp, r_a, above1) SELECT CASE (walltype) CASE (0) ! top cylinder call cyllweight(z,r,abovetmp, r_b, above2) CASE DEFAULT ! Ellipse as in gt170 call ellipseweight(z,r,abovetmp,r_0,z_0,invr_z,invr_r,Interior) END SELECT ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(:,0)=(Phidown*abovetmp(:,0) + Phiup*belowtmp(:,0) ) / & & (abovetmp(:,0)+belowtmp(:,0)) If(size(gtilde,2) .gt. 2) then ! first derivative denom=(abovetmp(:,0)+belowtmp(:,0))**2 gtilde(:,1)=(Phiup-Phidown)*(belowtmp(:,1)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,1)) / & & denom ! Axial derivative gtilde(:,2)=(Phiup-Phidown)*(belowtmp(:,2)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,2)) / & & denom ! Radial derivative End If End subroutine SUBROUTINE gUpDown(z,r,gtilde) ! standard g function added to rhs of poisson solver Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) Real(kind=db):: belowtmp(size(r,1),0:size(gtilde,2)-1), abovetmp(size(r,1),0:size(gtilde,2)-1) Real(kind=db):: denom(size(r,1)) ! Weight functions necessary for calculation of g call total_weight(z,r,belowtmp,abovetmp) ! Extension to the whole domain of the boundary conditions ! constructed by weight multiplied with boundary value gtilde(:,0)=(Phidown*abovetmp(:,0) + Phiup*belowtmp(:,0) ) / & & (abovetmp(:,0)+belowtmp(:,0)) If(size(gtilde,2) .gt. 2) then ! first derivative denom=(abovetmp(:,0)+belowtmp(:,0))**2 gtilde(:,1)=(Phiup-Phidown)*(belowtmp(:,1)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,1)) / & & denom ! Axial derivative gtilde(:,2)=(Phiup-Phidown)*(belowtmp(:,2)*abovetmp(:,0)-belowtmp(:,0)*abovetmp(:,2)) / & & denom ! Radial derivative End If END SUBROUTINE gUpDown Subroutine gtest(z,r,gtilde) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: gtilde(:,0:) ! 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(:), rsltarray(:) Integer:: fullrank, reducedrank Integer:: i,j call to_mat(full) fullrank=full%rank reducedrank=nbreducedspline call init(reducedrank,2,reduced) call init(fullrank,2,tempmat) !$OMP parallel private(j), firstprivate(temparray,rsltarray), default(shared) allocate(temparray(fullrank),rsltarray(fullrank)) !$OMP do do i=1,fullrank call getcol(etildet,i,temparray) rsltarray=vmx_mumps_mat_loc(full,temparray) DO j=1,fullrank CALL putele_sploc(tempmat%mat%row(i), j, rsltarray(j), .false.) END DO end do !$OMP end do !$OMP barrier !$OMP single call to_mat(tempmat) !$OMP end single !$OMP do do i=1,reducedrank call getcol(etildet,i,temparray) rsltarray=vmx_mumps_mat_loc(tempmat,temparray) DO j=1,reducedrank CALL putele_sploc(reduced%mat%row(i), j, rsltarray(j), .false.) END DO end do !$OMP end do deallocate(temparray,rsltarray) !$OMP end parallel end subroutine !!############################################################################################### ! subroutines imported from mumps_bsplines and sparse lib to have thread-private routines SUBROUTINE putele_sploc(arow, j, val, nlforce_zero) ! ! Put (overwrite) element j of row arow or insert it in an increasing "index" ! USE sparse TYPE(sprow), TARGET :: arow INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(in) :: val LOGICAL, INTENT(in), OPTIONAL :: nlforce_zero ! TYPE(elt), TARGET :: pre_root TYPE(elt), POINTER :: t, p LOGICAL :: rmnode ! pre_root%next => arow%row0 ! pre_root is linked to the head of the list. t => pre_root ! ! Remove node which has zero val or not? ! But never create new node with zero val ! rmnode = .TRUE. IF(PRESENT(nlforce_zero)) rmnode = .NOT.nlforce_zero ! DO WHILE( ASSOCIATED(t%next) ) p => t%next IF( p%index .EQ. j ) THEN IF(ABS(val).LE.EPSILON(0.0d0) .AND. rmnode) THEN ! Remove the node for zero val! t%next => p%next arow%nnz = arow%nnz-1 arow%row0 => pre_root%next ! In case the head is altered DEALLOCATE(p) ELSE p%val = val END IF RETURN END IF IF( p%index .GT. j ) EXIT t => t%next END DO ! ! Never create new node with zero val ! IF(ABS(val).GT.EPSILON(0.0d0)) THEN ALLOCATE(p) p = elt(j, val, t%next) t%next => p arow%nnz = arow%nnz+1 arow%row0 => pre_root%next END IF END SUBROUTINE putele_sploc ! SUBROUTINE getcol_loc(amat, j, arr) ! ! Get a column from sparse matrix ! USE mumps_bsplines TYPE(mumps_mat), INTENT(in) :: amat INTEGER, INTENT(in) :: j DOUBLE PRECISION, INTENT(out) :: arr(:) INTEGER :: i ! DO i=amat%istart,amat%iend CALL getele_loc(amat, i, j, arr(i)) END DO END SUBROUTINE getcol_loc SUBROUTINE getele_loc(mat, i, j, val) ! ! Get element (i,j) of sparse matrix ! USE mumps_bsplines TYPE(mumps_mat) :: mat INTEGER, INTENT(in) :: i, j DOUBLE PRECISION, INTENT(out) :: val INTEGER :: iget, jget INTEGER :: k, s, e ! iget = i jget = j IF(mat%nlsym) THEN IF( i.GT.j ) THEN ! Lower triangular part iget = j jget = i END IF END IF ! val = 0.0d0 ! Assume zero val if outside IF( iget.LT.mat%istart .OR. iget.GT.mat%iend ) RETURN ! IF(ASSOCIATED(mat%mat)) THEN CALL getele(mat%mat, iget, jget, val) ELSE s = mat%irow(iget) - mat%nnz_start + 1 e = mat%irow(iget+1) - mat%nnz_start k = isearch(mat%cols(s:e), jget) IF( k.GE.0 ) THEN val =mat%val(s+k) ELSE val = 0.0d0 ! Assume zero val if not found END IF END IF END SUBROUTINE getele_loc !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION vmx_mumps_mat_loc(mat, xarr) RESULT(yarr) ! ! Return product mat*x ! Use mumps_bsplines TYPE(mumps_mat) :: mat DOUBLE PRECISION, INTENT(in) :: xarr(:) DOUBLE PRECISION :: yarr(SIZE(xarr)) DOUBLE PRECISION :: alpha=1.0d0, beta=0.0d0 CHARACTER(len=6) :: matdescra INTEGER :: n ! n = mat%rank ! !#ifdef MKL IF(mat%nlsym) THEN matdescra = 'sun' ELSE matdescra = 'g' END IF ! CALL mkl_dcsrmv('N', n, n, alpha, matdescra, mat%val, & & mat%cols, mat%irow(1), mat%irow(2), xarr, & & beta, yarr) !#else ! yarr = 0.0d0 ! DO i=1,n ! DO j=mat%irow(i), mat%irow(i+1)-1 ! yarr(i) = yarr(i) + mat%val(j)*xarr(mat%cols(j)) ! END DO ! IF(mat%nlsym) THEN ! DO j=mat%irow(i)+1, mat%irow(i+1)-1 ! yarr(mat%cols(j)) = yarr(mat%cols(j)) & ! & + mat%val(j)*xarr(i) ! END DO ! END IF ! END DO ! Write(*,*) "linear mult") !#endif ! END FUNCTION vmx_mumps_mat_loc + END MODULE geometry \ No newline at end of file