diff --git a/src/weighttypes_mod.f90 b/src/weighttypes_mod.f90 index 5a2d6fe..414eb07 100644 --- a/src/weighttypes_mod.f90 +++ b/src/weighttypes_mod.f90 @@ -1,525 +1,526 @@ MODULE weighttypes USE constants use splinebound IMPLICIT NONE 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,z_a=0,z_b=0 Real(kind=db),save:: r_bLeft=0, r_bRight=0 Real(kind=db), save:: Phidown=0,Phiup=0 Integer, save::Interior=-1 Integer, save:: above1=1, above2=-1 type(spline_domain):: the_domain contains SUBROUTINE geom_w2(z,r,w,wupper) Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: walltmp(size(w,1),size(w,2)), elliptmp(size(w,1),size(w,2)) 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(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: belowtmp(size(w,1),size(w,2)), abovetmp(size(w,1),size(w,2)) ! 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(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: belowtmp(size(w,1),size(w,2)), abovetmp(size(w,1),size(w,2)) ! 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(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)) ! 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 ! 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(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! left metallic boundary call discweight(z,r, w1, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) 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(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! Right Dirichlet call discweight(z,r, w1, zgrid(nz),-1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! Left Dirichlet call discweight(z,r, w1, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w7 SUBROUTINE geom_w8(z,r,w, wupper) use basic, ONLY: zgrid, nz ! Defines the geometric weight for a central cylinder and a tilted upper cylinder with elliptic region ! and left and right dirichlet boundaries Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)) ! tilted wall call tiltedwall(z,r,w2,r_0,z_0, alpha, above2) ! calculate upper right boundary call cyllweight(z,r,w1,max(r_bRight,r_bLeft),above2) ! Intersection between slanted wall and upper right wall call Combine(w1,w2, w3, -1) ! calculate upper left boundary call cyllweight(z,r,w1,min(r_bRight,r_bLeft),above2) ! Union between previous weight w3 and upper left wall call Combine(w1,w3, w2, 1) ! tilted ellipse call tiltedellipse(z,r,w1,r_0,z_0, invr_z, invr_r, Interior, alpha) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, Interior) ! Left Dirichlet call discweight(z,r, w1, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w3, w2, -1) ! calculate coaxial insert call cyllweight(z,r,w1,r_a,above1) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, -1) ! Right Dirichlet call discweight(z,r, w1, zgrid(nz),-1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w8 SUBROUTINE geom_w10(z,r,w, wupper) ! Defines the geometric weight for a coaxial configuration with end-plugs Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)), w4(size(w,1),size(w,2)) ! Right Dirichlet call discweight(z,r, w1, z_b,-1) ! Left Dirichlet call discweight(z,r, w2, z_a,1) ! Union between previous weight w3 and upper left wall call Combine(w1,w2, w3, -1) ! calculate upper radial boundary call cyllweight(z,r,w1,r_b,-1) ! calculate lower radial boundary call cyllweight(z,r,w2,r_a,1) ! Union between previous weight w3 and upper left wall call Combine(w1,w2, w4, -1) if( present(wupper) ) then wupper=w4 w=w3 else ! gives total weight call Combine(w3, w4, w, -1) end if End SUBROUTINE geom_w10 SUBROUTINE geom_w11(z,r,w, wupper) ! Defines the geometric weight for the inside of a ring Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(0:size(w,1)-1,size(w,2)), w2(0:size(w,1)-1,size(w,2)) ! elllipse centered in r_0, z_0 call ellipseweight2(z,r,w1,r_0,z_0, invr_z, invr_r, Interior) ! elllipse centered in r_0, z_0+z_r/2 call ellipseweight2(z,r,w2,r_0,z_0+z_r/3, 3*invr_z, 3*invr_r, -Interior) if( present(wupper) ) then wupper=w2 w=w1 else ! gives total weight call Combine(w2, w1, w, -1) end if End SUBROUTINE geom_w11 SUBROUTINE geom_w12(z,r,w, wupper) Use basic, only: ZGRID ! Defines the geometric weight for a central cylinder with elliptic cut ! and left and right dirichlet boundaries Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) Real(kind=db):: w1(size(w,1),size(w,2)), w2(size(w,1),size(w,2)), w3(size(w,1),size(w,2)) Real(kind=db):: w4(size(w,1),size(w,2)) ! calculate outer cylinder call cyllweight(z,r,w1,r_b,above2) ! Left Dirichlet call discweight(z,r, w2, zgrid(0),1) ! Intersection of ellipse and previous weight call Combine(w1, w2, w3, -1) ! calculate coaxial insert call cyllweight(z,r,w2,r_a,above1) ! calculate ellipse cut on the inside call ellipseweight2(z,r,w4, r_a, z_0, invr_z, invr_r, 1) call Combine(w2,w4,w1,1) if( present(wupper) ) then wupper=w3 w=w1 else ! gives total weight call Combine(w1, w3, w, -1) end if End SUBROUTINE geom_w12 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,1) .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,1) .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,1) .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_c, z_c, invrz, invrr, Inside) Real(kind=db):: z(:),r(:),w(:,:) Integer:: Inside, above Real(kind=db):: r_lim, r_c,z_c,invrz,invrr Real(kind=db):: walltmp(size(w,1),size(w,2)), elliptmp(size(w,1),size(w,2)) call cyllweight(z,r,walltmp, r_lim, above) call ellipseweight2(z,r,elliptmp,r_c, z_c, invrz, invrr, Inside) call combine( walltmp, elliptmp, w, Inside) END SUBROUTINE ellipsewall Subroutine ellipseweight(z,r,w, r_c, z_c, invrz, invrr, Inside) Real(kind=db):: z(:),r(:),w(0:,:) Real(kind=db):: r_c,z_c,invrz,invrr Integer:: Inside Real(kind=db):: D(size(r,1)) D=sqrt((r-r_c)**2*invrr**2+(z-z_c)**2*invrz**2) w(0,:)=Inside*(1-D) ! weight at position r,z If(size(w,1) .gt. 1) then ! first derivative w(1,:)=-(z-z_c)*invrz**2/D*Inside ! z derivative of w w(2,:)=-(r-r_c)*invrr**2/D*Inside ! r derivative of w End If End subroutine ellipseweight Subroutine ellipseweight2(z,r,w, r_c, z_c, invrz, invrr, Inside) Real(kind=db):: z(:),r(:),w(0:,:) Real(kind=db):: r_c,z_c,invrz,invrr Integer:: Inside Real(kind=db):: D(size(r,1)) D=(r-r_c)**2*invrr**2+(z-z_c)**2*invrz**2 w(0,:)=Inside*(1-D) ! weight at position r,z If(size(w,1) .gt. 1) then ! first derivative w(1,:)=-2*(z-z_c)*invrz**2*Inside ! z derivative of w w(2,:)=-2*(r-r_c)*invrr**2*Inside ! r derivative of w End If + w=w/max(invrz,invrr) End subroutine ellipseweight2 Subroutine tiltedellipse(z,r,w, r_c, z_c, invrz, invrr, Inside, alpha) Real(kind=db):: z(:),r(:),w(0:,:) Real(kind=db):: r_c,z_c,invrz,invrr, alpha, cosa, sina Real(kind=db):: deltar(size(r,1)), deltaz(size(z,1)) Integer:: Inside Real(kind=db):: D(size(r,1)) cosa=cos(alpha) sina=sin(alpha) deltar=(r-r_c) deltaz=(z-z_c) D=sqrt(((deltaz*cosa+deltar*sina)*invrz)**2+((deltaz*sina-deltar*cosa)*invrr)**2) w(0,:)=1-D ! weight at position r,z If(size(w,1) .gt. 1) then ! first derivative w(1,:)=(-cosa*(deltaz*cosa+deltar*sina)*invrz**2-sina*(deltaz*sina-deltar*cosa)*invrr**2)/D ! z derivative of w w(2,:)=(-sina*(deltaz*cosa+deltar*sina)*invrz**2+cosa*(deltaz*sina-deltar*cosa)*invrr**2)/D ! r derivative of w End If w=w*Inside 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,2)), denom(size(w,2)) 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,1) .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 geom_spline(z,r,w,wupper) Use splinebound, ONLY: spline_w Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) Real(kind=db), OPTIONAL:: wupper(:,:) call spline_w(the_domain,z,r,w) End SUBROUTINE geom_spline SUBROUTINE geom_splinetot(z,r,w,idwall) Use splinebound, ONLY: spline_w Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: w(:,:) INTEGER, optional, INTENT(OUT):: idwall(:) call spline_wtot(the_domain,z,r,w,idwall) End SUBROUTINE geom_splinetot SUBROUTINE gspline(z,r,g,w) Use splinebound, ONLY: spline_g Real(kind=db), INTENT(IN):: r(:),z(:) Real(kind=db), INTENT(OUT):: g(:,:) Real(kind=db), INTENT(IN),OPTIONAL::w(:,:) call spline_g(the_domain,z,r,g,w) End SUBROUTINE gspline END MODULE weighttypes