-
Notifications
You must be signed in to change notification settings - Fork 37
Expand file tree
/
Copy pathdiffusion.F90
More file actions
170 lines (127 loc) · 4.45 KB
/
Copy pathdiffusion.F90
File metadata and controls
170 lines (127 loc) · 4.45 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
module diffusion
!! This module sets up the ionospheric diffusion problem and then passes it off to the parabolic solvers.
use phys_consts, only: kB,lsp,gammas,mindensdiv, wp
!! sizes are not imported in case we want to do subgrid diffusion
use mesh, only : curvmesh
use grid, only : gridflag
use PDEparabolic, only : backEuler1D,TRBDF21D
implicit none
interface TRBDF23D
module procedure TRBDF23D_curv
end interface TRBDF23D
interface backEuler3D
module procedure backEuler3D_curv
end interface backEuler3D
private
public :: diffusion_prep, TRBDF23D, backEuler3D
contains
pure subroutine diffusion_prep(isp,x,lambda,betacoeff,ns,T,A,B,C,D,E,Tn,Teinf)
!! COMPUTE COEFFICIENTS IN DIFFUSION EQUATIONS AND LOAD UP GHOST CELLS
!!
!! Note: done on a per species basis. This is never called over the full grid
integer, intent(in) :: isp
type(curvmesh), intent(in) :: x
real(wp), dimension(:,:,:), intent(in) :: lambda,betacoeff
real(wp), dimension(-1:,-1:,-1:), intent(in) :: ns
real(wp), dimension(-1:,-1:,-1:), intent(inout) :: T
real(wp), dimension(1:size(ns,1)-4,1:size(ns,2)-4,1:size(ns,3)-4), &
intent(out) :: A,B,C,D,E
real(wp), dimension(:,:,:), intent(in) :: Tn
real(wp), intent(in) :: Teinf
real(wp) :: Tn0
integer :: lx1,lx2,lx3,ix2,ix3
lx1=size(ns,1)-4
lx2=size(ns,2)-4
lx3=size(ns,3)-4
!COEFFICIENTS OF PARABOLIC EQUAITON
A(:,:,:)=0._wp
C(:,:,:)=(gammas(isp)-1._wp)/kB/max(ns(1:lx1,1:lx2,1:lx3),mindensdiv)/ &
(x%h1(1:lx1,1:lx2,1:lx3)*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3))
B(:,:,:)=C(:,:,:)*betacoeff/x%h1(1:lx1,1:lx2,1:lx3) !beta must be set to zero if not electrons!
D=lambda*x%h2(1:lx1,1:lx2,1:lx3)*x%h3(1:lx1,1:lx2,1:lx3)/x%h1(1:lx1,1:lx2,1:lx3)
E=0._wp
!SET THE BOUNDARY CONDITIONS BASED ON GRID TYPE
if (gridflag==0) then
do ix3=1,lx3
do ix2=1,lx2
Tn0=Tn(1,ix2,ix3)
T(0,ix2,ix3)=Tn0
Tn0=Tn(lx1,ix2,ix3)
T(lx1+1,ix2,ix3)=Tn0
end do
end do
else if (gridflag==1) then
do ix3=1,lx3
do ix2=1,lx2
Tn0=Tn(lx1,ix2,ix3)
T(lx1+1,ix2,ix3)=Tn0 !bottom
T(0,ix2,ix3)=Teinf !top
end do
end do
else
do ix3=1,lx3
do ix2=1,lx2
Tn0=Tn(1,ix2,ix3)
T(0,ix2,ix3)=Tn0 !bottom
T(lx1+1,ix2,ix3)=Teinf !top
end do
end do
end if
end subroutine diffusion_prep
function backEuler3D_curv(f,A,B,C,D,E,dt,x)
!! SOLVE A 3D SEQUENCE OF 1D DIFFUSION PROBLEMS.
!! GHOST CELLS ARE ACCOMODATED AS THEY PROVIDE
!! A CONVENIENT MEMORY SPACE FOR BOUNDARY CONDITIONS.
real(wp), dimension(:,:,:), intent(in) :: A,B,C,D,E !trimmed to grid size
real(wp), dimension(-1:,-1:,-1:), intent(in) :: f !expected to include ghosts
real(wp), intent(in) :: dt
type(curvmesh), intent(in) :: x
integer :: ix1,ix2,ix3,lx1,lx2,lx3
real(wp),dimension(size(f,1)-4) :: fx1slice
real(wp), dimension(-1:size(f,1)-2,-1:size(f,2)-2,-1:size(f,3)-2) :: backEuler3D_curv
lx1=size(f,1)-4
lx2=size(f,2)-4
lx3=size(f,3)-4
do ix3=1,lx3
do ix2=1,lx2
fx1slice=f(1:lx1,ix2,ix3)
fx1slice=backEuler1D(fx1slice,A(:,ix2,ix3), &
B(:,ix2,ix3),C(:,ix2,ix3),D(:,ix2,ix3),E(:,ix2,ix3), &
f(0,ix2,ix3),f(lx1+1,ix2,ix3),dt,x%dx1,x%dx1i)
!! inner ghost cells include boundary conditions
backEuler3D_curv(1:lx1,ix2,ix3)=fx1slice
end do
end do
end function backEuler3D_curv
function TRBDF23D_curv(f,A,B,C,D,E,dt,x)
!! SOLVE A 3D SEQUENCE OF 1D DIFFUSION PROBLEMS.
!! GHOST CELLS ARE ACCOMODATED AS THEY PROVIDE
!! A CONVENIENT MEMORY SPACE FOR BOUNDARY CONDITIONS.
!!
!! Note that this function also plays the role of abstracting
!! away the grid structure so that the individual 1D lines are
!! solved using arrays for differences instead of structure members
!> trimmed to grid size
real(wp), dimension(:,:,:), intent(in) :: A,B,C,D,E
!> expected to include ghosts
real(wp), dimension(-1:,-1:,-1:), intent(in) :: f
real(wp), intent(in) :: dt
type(curvmesh), intent(in) :: x
integer :: ix1,ix2,ix3,lx1,lx2,lx3
real(wp),dimension(size(f,1)-4) :: fx1slice
real(wp), dimension(-1:size(f,1)-2,-1:size(f,2)-2,-1:size(f,3)-2) :: TRBDF23D_curv
lx1=size(f,1)-4
lx2=size(f,2)-4
lx3=size(f,3)-4
do ix3=1,lx3
do ix2=1,lx2
fx1slice=f(1:lx1,ix2,ix3)
fx1slice=TRBDF21D(fx1slice,A(:,ix2,ix3), &
B(:,ix2,ix3),C(:,ix2,ix3),D(:,ix2,ix3),E(:,ix2,ix3), &
f(0,ix2,ix3),f(lx1+1,ix2,ix3),dt,x%dx1,x%dx1i)
!! inner ghost cells include boundary conditions
TRBDF23D_curv(1:lx1,ix2,ix3)=fx1slice
end do
end do
end function TRBDF23D_curv
end module diffusion