ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
equilibrium_augmenter.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 module equilibrium_augmenter
9 
10  use euitm_schemas
11  use itm_types
12  use itm_constants
13  use deallocate_structures
14 
15 contains
16 
17  subroutine augment_psi_rz(equilibrium)
18 
19  implicit none
20 
21  type (type_equilibrium) :: equilibrium
22  integer, parameter :: nr=65, nz=129
23  integer :: npsi, ndim1, ndim2
24  integer :: idim1, idim2
25 
26  real (R8), allocatable :: psi_in(:,:), f_dia_in(:,:), w(:,:)
27 
28  real (R8), allocatable :: r(:), z(:), ans(:), tx(:), ty(:), c(:), &
29  wrk1(:), wrk2(:), wrk(:)
30  integer, allocatable :: iwrk(:), jwrk(:)
31 
32  real (R8) :: rmin, rmax, zmin, zmax, dr, dz, fp, s, r0, b0, psi_ax, psi_bnd, sign_psi
33  integer :: kx, ky, m,nxest, nyest, nmax, u, v, km, ne, bx, by, &
34  b1, b2, lwrk1, lwrk2, lwrk, kwrk, mwrk, nx, ny, ier, i, j
35 
36 
37  npsi = size(equilibrium%profiles_1d%psi, dim=1)
38  ndim1 = size(equilibrium%coord_sys%position%R, dim=1)
39  ndim2 = size(equilibrium%coord_sys%position%R, dim=2)
40 
41  if(npsi.ne.ndim1) then
42  write(*,*) 'NPSI <> NDIM1: ',npsi, ndim1
43  stop 'NPSI <> NDIM1'
44  endif
45 
46  r0=equilibrium%global_param%toroid_field%r0
47  b0=equilibrium%global_param%toroid_field%b0
48  psi_ax=equilibrium%profiles_1d%psi(1)
49  psi_bnd=equilibrium%profiles_1d%psi(npsi)
50  sign_psi=sign(1.0_r8, psi_bnd-psi_ax)
51 
52  allocate(psi_in(ndim1,ndim2), f_dia_in(ndim1,ndim2), w(ndim1,ndim2))
53  allocate(r(nr), z(nz))
54 
55  do idim2=1, ndim2
56  psi_in(:,idim2) = equilibrium%profiles_1d%psi(:)
57  f_dia_in(:,idim2) = equilibrium%profiles_1d%f_dia(:)
58  enddo
59 
60  w=1.0_r8
61 
62  rmin=minval(equilibrium%coord_sys%position%R)
63  rmax=maxval(equilibrium%coord_sys%position%R)
64  zmin=minval(equilibrium%coord_sys%position%Z)
65  zmax=maxval(equilibrium%coord_sys%position%Z)
66 
67  dr=rmax-rmin
68  dz=zmax-zmin
69 
70  rmin=rmin-dr/20.0_r8
71  rmax=rmax+dr/20.0_r8
72  zmin=zmin-dz/20.0_r8
73  zmax=zmax+dz/20.0_r8
74 
75  dr=(rmax-rmin)/(nr-1)
76  dz=(zmax-zmin)/(nz-1)
77 
78  do i=1, nr
79  r(i)=rmin+(i-1)*dr
80  enddo
81  do i=1, nz
82  z(i)=zmin+(i-1)*dz
83  enddo
84 
85  kx=3
86  ky=3
87 
88  m = ndim1*ndim2
89  nxest = max(2*(kx+1),kx+1+int(sqrt(m/2.0_r8)))
90  nyest = max(2*(ky+1),ky+1+int(sqrt(m/2.0_r8)))
91  nmax = max(nxest, nyest)
92  u = nxest-kx-1
93  v = nyest-ky-1
94  km = max(kx,ky)+1
95  ne = max(nxest,nyest)
96  bx = kx*v+ky+1
97  by = ky*u+kx+1
98  if(bx.le.by) then
99  b1 = bx
100  b2 = b1+v-ky
101  else
102  b1 = by
103  b2 = b1+u-kx
104  endif
105  lwrk1 = u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
106  lwrk2 = u*v*(b2+1)+b2
107  kwrk = m+(nxest-2*kx-1)*(nyest-2*ky-1)
108 
109  allocate(tx(nmax), ty(nmax), c((nxest-kx-1)*(nyest-ky-1)), &
110  wrk1(lwrk1), wrk2(lwrk2), iwrk(kwrk))
111 
112  s = sqrt(real(m,r8))/100
113 
114 ! spline psi
115  call surfit(0, m, &
116  equilibrium%coord_sys%position%R, equilibrium%coord_sys%position%Z, psi_in, &
117  w, rmin, rmax, zmin, zmax, kx, ky, s, nxest, nyest, nmax, 1e-15_r8, &
118  nx, tx, ny, ty, c, fp, wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
119 
120  write(*,*) 'surfit returned ier = ', ier
121  write(*,*) 'nx, ny, fp = ', nx, ny, fp
122 
123  lwrk = nr*(kx+1)+nz*(ky+1)+(nx-1)*(ny-1)
124  mwrk = nr+nz
125 
126  allocate(wrk(lwrk), jwrk(mwrk), ans(nr*nz))
127 
128  call bispev(tx, nx, ty, ny, c, kx, ky, r, nr, z, nz, ans, wrk, lwrk, jwrk, mwrk, ier)
129 
130  write(*,*) 'bispev returned ier = ', ier
131 
132  if(associated(equilibrium%profiles_2d)) call deallocate_cpo(equilibrium%profiles_2d)
133  allocate(equilibrium%profiles_2d(1))
134 
135  allocate(equilibrium%profiles_2d(1)%psi(nr,nz))
136  do i=1, nr
137  do j=1, nz
138  equilibrium%profiles_2d(1)%psi(i,j) = ans(nz*(i-1)+j)
139  enddo
140  enddo
141 
142 ! br
143  call parder(tx, nx, ty, ny, c, kx, ky, 0, 1, r, nr, z, nz, ans, wrk, lwrk, jwrk, mwrk, ier)
144 
145  write(*,*) 'parder returned ier = ', ier
146 
147  allocate(equilibrium%profiles_2d(1)%br(nr,nz))
148  do i=1, nr
149  do j=1, nz
150  equilibrium%profiles_2d(1)%br(i,j) = -ans(nz*(i-1)+j)/r(i)/2.0_r8/itm_pi
151  enddo
152  enddo
153 
154 !bz
155  call parder(tx, nx, ty, ny, c, kx, ky, 1, 0, r, nr, z, nz, ans, wrk, lwrk, jwrk, mwrk, ier)
156 
157  write(*,*) 'parder returned ier = ', ier
158 
159  allocate(equilibrium%profiles_2d(1)%bz(nr,nz))
160  do i=1, nr
161  do j=1, nz
162  equilibrium%profiles_2d(1)%bz(i,j) = ans(nz*(i-1)+j)/r(i)/2.0_r8/itm_pi
163  enddo
164  enddo
165 
166  s = sqrt(real(m,r8))/10000
167 
168 ! spline f_dia
169  call surfit(0, m, &
170  equilibrium%coord_sys%position%R, equilibrium%coord_sys%position%Z, f_dia_in, &
171  w, rmin, rmax, zmin, zmax, kx, ky, s, nxest, nyest, nmax, 1e-15_r8, &
172  nx, tx, ny, ty, c, fp, wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
173 
174  write(*,*) 'surfit returned ier = ', ier
175  write(*,*) 'nx, ny, fp = ', nx, ny, fp
176 
177  deallocate(wrk, jwrk, ans)
178 
179  lwrk = nr*(kx+1)+nz*(ky+1)+(nx-1)*(ny-1)
180  mwrk = nr+nz
181 
182  allocate(wrk(lwrk), jwrk(mwrk), ans(nr*nz))
183 
184  call bispev(tx, nx, ty, ny, c, kx, ky, r, nr, z, nz, ans, wrk, lwrk, jwrk, mwrk, ier)
185 
186  write(*,*) 'bispev returned ier = ', ier
187 
188  allocate(equilibrium%profiles_2d(1)%bphi(nr,nz))
189  do i=1, nr
190  do j=1, nz
191  if(equilibrium%profiles_2d(1)%psi(i,j)*sign_psi .gt. psi_bnd*sign_psi) then
192  equilibrium%profiles_2d(1)%bphi(i,j) = r0*b0/r(i)
193  else
194  equilibrium%profiles_2d(1)%bphi(i,j) = ans(nz*(i-1)+j)/r(i)
195  endif
196  enddo
197  enddo
198 
199  allocate(equilibrium%profiles_2d(1)%grid_type(4))
200  allocate(equilibrium%profiles_2d(1)%grid%dim1(nr))
201  allocate(equilibrium%profiles_2d(1)%grid%dim2(nz))
202  equilibrium%profiles_2d(1)%grid_type(1)='1'
203  equilibrium%profiles_2d(1)%grid_type(2)='rectangular'
204  equilibrium%profiles_2d(1)%grid_type(3)='0'
205  equilibrium%profiles_2d(1)%grid_type(4)='not relevant'
206  equilibrium%profiles_2d(1)%grid%dim1=r
207  equilibrium%profiles_2d(1)%grid%dim2=z
208 
209  return
210 
211  end subroutine augment_psi_rz
212 
213 end module equilibrium_augmenter
bispev
subroutine bispev(tx, nx, ty, ny, c, kx, ky, x, mx, y, my, z, wrk, lwrk, iwrk, kwrk, ier)
Definition: bispev.f:1
equilibrium_augmenter::augment_psi_rz
subroutine augment_psi_rz(equilibrium)
Definition: equilibrium_augmenter.f90:17
phys_functions::b2
real(r8) function b2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:289
equilibrium_augmenter
Augment an inverse equilibrium with psi, Br, Bz and Bphi as a function of R and Z.
Definition: equilibrium_augmenter.f90:8
surfit
subroutine surfit(iopt, m, x, y, z, w, xb, xe, yb, ye, kx, ky, s, nxest, nyest, nmax, eps, nx, tx, ny, ty, c, fp, wrk1, lwrk1, wrk2, lwrk2, iwrk, kwrk, ier)
Definition: surfit.f:1
parder
subroutine parder(tx, nx, ty, ny, c, kx, ky, nux, nuy, x, mx, y, my, z, wrk, lwrk, iwrk, kwrk, ier)
Definition: parder.f:1