Homec4science
Diffusion FENNECS: Finite Element Non-Neutral Electron Cloud Simulator (geom_test) Browse

fennecs/src/weighttypes_mod.f90
c687afcc43d6geom_test

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
Real(kind=db), save:: L_r=0.16, L_z=0.24
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(:,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=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(:,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)
! 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)
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), w4(size(w,1),0:size(w,2)-1)
! 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)
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)
! elllipse centered in r_0, z_0
call ellipseweight2(z,r,w1,r_0,z_0, invr_z, invr_r, Interior)
w=w1
if( present(wupper) ) then
wupper=0
end if
End SUBROUTINE geom_w11
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_c, z_c, invrz, invrr, Inside)
Real(kind=db):: z(:),r(:),w(:,0:)
Integer:: Inside, above
Real(kind=db):: r_lim, r_c,z_c,invrz,invrr
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_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,2) .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,2) .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
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,2) .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,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 geom_spline(z,r,w,wupper)
Use splinebound, ONLY: spline_w
Real(kind=db), INTENT(IN):: r(:),z(:)
Real(kind=db), INTENT(OUT):: w(:,0:)
Real(kind=db), OPTIONAL:: wupper(:,0:)
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(:,0:)
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(:,0:)
Real(kind=db), INTENT(IN),OPTIONAL::w(:,0:)
call spline_g(the_domain,z,r,g,w)
End SUBROUTINE gspline
END MODULE weighttypes
c4science · Help