Homec4science
Diffusion FENNECS lite (master) Browse

fennecslite/src/weighttypes_mod.f90
b8eece80d5a5master

Last Change
Raw File

weighttypes_mod.f90
View Options

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(r,z,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
c4science · Help