ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
wrapper.F90
Go to the documentation of this file.
1 Module auxiliaries
2 
3 CONTAINS
4 
5 SUBROUTINE prof_init(eq, coreprof)
6 
7  USE phys_constants
8  USE euitm_schemas
9 
10  IMPLICIT NONE
11 
12  TYPE (type_equilibrium) :: eq
13  TYPE (type_coreprof) :: coreprof
14 
15  INTEGER(ITM_I4) :: i
16  INTEGER(ITM_I4) :: nrho_prof = 20, nion = 1, npsi_eq = 33, neta_eq = 64
17 
18  REAL(R8), DIMENSION(:), POINTER :: eta,ra,psi,qq,jphi,nne,tte
19 
20  REAL(R8) :: r_minor=0.5_r8, r_major=1.65_r8, mag_field=-2.5_r8, &
21  q_axis = -1.5_r8, j_axis, psi_axis
22 
23 !... load profiles
24 
25  ALLOCATE(coreprof%composition%amn(nion))
26  ALLOCATE(coreprof%composition%zion(nion))
27 
28  ALLOCATE(coreprof%rho_tor_norm(nrho_prof))
29  ALLOCATE(coreprof%rho_tor(nrho_prof))
30  ALLOCATE(coreprof%psi%value(nrho_prof))
31  ALLOCATE(coreprof%drho_dt(nrho_prof))
32 
33  ALLOCATE(coreprof%ne%value(nrho_prof))
34  ALLOCATE(coreprof%te%value(nrho_prof))
35  ALLOCATE(coreprof%ni%value(nrho_prof,nion))
36  ALLOCATE(coreprof%ti%value(nrho_prof,nion))
37 
38  ALLOCATE(coreprof%profiles1d%q%value(nrho_prof))
39  ALLOCATE(coreprof%profiles1d%zeff%value(nrho_prof))
40  ALLOCATE(coreprof%profiles1d%wtor%value(nrho_prof,nion))
41  ALLOCATE(coreprof%profiles1d%pe%value(nrho_prof))
42  ALLOCATE(coreprof%profiles1d%pi%value(nrho_prof,nion))
43  ALLOCATE(coreprof%profiles1d%pr_th%value(nrho_prof))
44  ALLOCATE(coreprof%profiles1d%jtot%value(nrho_prof))
45  ALLOCATE(coreprof%profiles1d%jni%value(nrho_prof))
46  ALLOCATE(coreprof%profiles1d%joh%value(nrho_prof))
47  ALLOCATE(coreprof%profiles1d%qoh%value(nrho_prof))
48  ALLOCATE(coreprof%profiles1d%sigmapar%value(nrho_prof))
49 
50 !... set the global parameters
51 
52  coreprof%toroid_field%b0 = mag_field
53  coreprof%toroid_field%r0 = r_major
54 
55  coreprof%composition%amn = 1.0_r8
56  coreprof%composition%zion = 1.0_r8
57 
58  eq%global_param%toroid_field%b0 = mag_field
59  eq%global_param%toroid_field%r0 = r_major
60 
61  eq%global_param%mag_axis%q = q_axis
62 
63  eq%eqgeometry%geom_axis%r = r_major
64  eq%eqgeometry%geom_axis%z = 0._r8
65  eq%eqgeometry%a_minor = r_minor
66  eq%eqgeometry%elongation = 1.0_r8
67  eq%eqgeometry%tria_upper = 0._r8
68  eq%eqgeometry%tria_lower = 0._r8
69 
70 !... load equil
71 
72  ALLOCATE(eq%codeparam%codename(1))
73  eq%codeparam%codename(1) = "wrapper input for BDSEQ"
74 
75  ALLOCATE(eq%profiles_1d%psi(npsi_eq))
76  ALLOCATE(eq%profiles_1d%phi(npsi_eq))
77  ALLOCATE(eq%profiles_1d%rho_tor(npsi_eq))
78  ALLOCATE(eq%profiles_1d%rho_vol(npsi_eq))
79  ALLOCATE(eq%profiles_1d%pressure(npsi_eq))
80  ALLOCATE(eq%profiles_1d%jparallel(npsi_eq))
81  ALLOCATE(eq%profiles_1d%jphi(npsi_eq))
82  ALLOCATE(eq%profiles_1d%pprime(npsi_eq))
83  ALLOCATE(eq%profiles_1d%ffprime(npsi_eq))
84  ALLOCATE(eq%profiles_1d%F_dia(npsi_eq))
85  ALLOCATE(eq%profiles_1d%q(npsi_eq))
86 
87  eq%profiles_1d%phi=0._r8
88  eq%profiles_1d%psi=0._r8
89  eq%profiles_1d%rho_tor=0._r8
90  eq%profiles_1d%rho_vol=0._r8
91  eq%profiles_1d%jphi=0._r8
92  eq%profiles_1d%jparallel=0._r8
93  eq%profiles_1d%pressure=0._r8
94  eq%profiles_1d%q=0._r8
95  eq%profiles_1d%pprime=0._r8
96  eq%profiles_1d%ffprime=0._r8
97  eq%profiles_1d%F_dia=0._r8
98 
99 !... eq grid
100 
101  ALLOCATE(eq%coord_sys%grid_type(1))
102  ALLOCATE(eq%coord_sys%grid%dim1(npsi_eq))
103  ALLOCATE(eq%coord_sys%grid%dim2(neta_eq))
104  eq%coord_sys%grid_type = "unit rho theta cyl grid for input"
105  eq%coord_sys%grid%dim1 = &
106  (/ ((1.0_r8/(npsi_eq-1))*REAL(i-1),i=1,npsi_eq) /)
107  eq%coord_sys%grid%dim2 = &
108  (/ ((tpi/neta_eq)*REAL(i-1),i=1,neta_eq) /)
109 
110  ra => eq%coord_sys%grid%dim1
111  eta => eq%coord_sys%grid%dim2
112 
113 !... plasma boundary
114 
115  ALLOCATE(eq%eqgeometry%boundary(1))
116  ALLOCATE(eq%eqgeometry%boundary(1)%r(neta_eq))
117  ALLOCATE(eq%eqgeometry%boundary(1)%z(neta_eq))
118 
119  eq%eqgeometry%boundary(1)%r = &
120  eq%eqgeometry%geom_axis%r + &
121  eq%eqgeometry%a_minor*cos(eta)
122  eq%eqgeometry%boundary(1)%z = &
123  eq%eqgeometry%geom_axis%z + &
124  eq%eqgeometry%a_minor*sin(eta)
125 
126 !... eq profiles
127 
128  qq => eq%profiles_1d%q
129  psi => eq%profiles_1d%psi
130  jphi => eq%profiles_1d%jphi
131 
132  eq%profiles_1d%rho_tor=eq%eqgeometry%a_minor * ra
133  eq%profiles_1d%rho_vol = ra
134 
135  j_axis = - 2.0_r8*mag_field/(mu_0*r_major*q_axis)
136  psi_axis = (pi*r_minor*r_minor)*mag_field/q_axis
137 
138  qq=q_axis*exp(ra*ra)
139  jphi = j_axis*exp(-ra*ra)*(1.0_r8-ra*ra)
140  psi = psi_axis*(1.0_r8-exp(-ra*ra))
141 
142  eq%profiles_1d%jparallel=eq%profiles_1d%jphi
143 
144  eq%profiles_1d%F_dia=r_major*mag_field
145 
146  eq%profiles_1d%phi=pi*mag_field*eq%profiles_1d%rho_tor*eq%profiles_1d%rho_tor
147 
148 !... coreprof grid and profiles
149 
150  coreprof%rho_tor_norm=(/ ((1.0_r8/nrho_prof)*(i-0.5),i=1,nrho_prof) /)
151  coreprof%rho_tor=coreprof%rho_tor_norm * eq%eqgeometry%a_minor
152 
153  ra => coreprof%rho_tor_norm
154  nne => coreprof%ne%value
155  tte => coreprof%te%value
156 
157  coreprof%ne%value=5.0e19_r8*exp(-8.0_r8*ra*ra/7.0_r8)
158  coreprof%te%value=2.0e3_r8*exp(-8.0_r8*ra*ra/3.0_r8)
159  coreprof%ni%value(:,1)=coreprof%ne%value
160  coreprof%ti%value(:,1)=coreprof%te%value
161  coreprof%profiles1d%zeff%value = 2.0_r8
162  coreprof%profiles1d%wtor%value=0.0_r8
163  coreprof%time=0.0_r8
164 
165  coreprof%profiles1d%pe%value=coreprof%ne%value*kb*coreprof%te%value
166  coreprof%profiles1d%pi%value=coreprof%ni%value*kb*coreprof%ti%value
167  coreprof%profiles1d%pr_th%value=coreprof%profiles1d%pe%value &
168  + coreprof%profiles1d%pi%value(:,1)
169 
170 !... put q and psi and J onto coreprof grid
171 
172  CALL l3interp( eq%profiles_1d%q, eq%profiles_1d%rho_tor, npsi_eq, &
173  coreprof%profiles1d%q%value, coreprof%rho_tor, nrho_prof)
174  CALL l3interp( eq%profiles_1d%psi, eq%profiles_1d%rho_tor, npsi_eq, &
175  coreprof%psi%value, coreprof%rho_tor, nrho_prof)
176  CALL l3interp( eq%profiles_1d%jphi, eq%profiles_1d%rho_tor, npsi_eq, &
177  coreprof%profiles1d%jtot%value, coreprof%rho_tor, nrho_prof)
178 
179  coreprof%profiles1d%jni%value=0.0_r8
180  coreprof%profiles1d%joh%value=coreprof%profiles1d%jtot%value
181 
182 !... put p onto equil grid and get derivs
183 
184  CALL l3interp( coreprof%profiles1d%pr_th%value, &
185  coreprof%rho_tor, nrho_prof, &
186  eq%profiles_1d%pressure, eq%profiles_1d%rho_tor, npsi_eq)
187 
188  CALL l3deriv( eq%profiles_1d%pressure, eq%profiles_1d%psi, npsi_eq, &
189  eq%profiles_1d%pprime, eq%profiles_1d%psi, npsi_eq)
190 
191  eq%profiles_1d%ffprime = mu_0*r_major*(jphi/tpi - r_major* &
192  eq%profiles_1d%pprime)
193 
194 !... set loop voltage
195 
196  coreprof%globalparam%vloop = 0._r8
197 
198 !... stamp the time
199 
200  coreprof%time=0.0_r8
201 
202  eq%time=0.0_r8
203 
204 END SUBROUTINE prof_init
205 END Module auxiliaries
206 
207 
208 PROGRAM bdseq_wrapper
209 
210  !... uses either simple input or a CPO file called IMP4Init to wrap BDSEQ
211  !... simple input if the file is absent
212 
213 #ifdef MPI
214  USE mpi
215 #endif
216 
217  USE itm_constants
218  USE euitm_schemas
219  USE read_structures
220 
221  USE auxiliaries
222 
223  IMPLICIT NONE
224 
225  INTEGER :: ios
226 
227  TYPE (type_equilibrium), pointer :: eq_in(:), eq(:)
228  TYPE (type_coreprof), POINTER :: coreprof(:)
229  TYPE (type_param) :: code_parameters
230 
231  interface
232  subroutine bdseq(eq_in, eq, code_parameters)
233  use euitm_schemas
234  type (type_equilibrium), pointer :: eq_in(:), eq(:)
235  type (type_param) :: code_parameters
236  end subroutine bdseq
237  end interface
238 
239 !... this gets a copy of the code onto every core
240 
241 #ifdef MPI
242  INTEGER :: ierr
243  CALL mpi_init(ierr)
244 #endif
245 
246 !... set initial state
247 
248  ALLOCATE (eq_in(1))
249  ALLOCATE (coreprof(1))
250 
251  OPEN (unit = 10, file = 'IMP4Init', &
252  status = 'old', form = 'formatted', &
253  action = 'read', iostat = ios)
254 
255  IF (ios == 0) THEN
256 
257  CLOSE (10)
258 
259  call open_read_file(10, 'IMP4Init' )
260  call read_cpo(coreprof(1), 'IMP4InitCoreprof' )
261  call read_cpo(eq_in(1), 'IMP4InitEquil' )
262  call close_read_file
263 
264  ELSE
265 
266  call prof_init(eq_in(1),coreprof(1))
267 
268  END IF
269 
270 !... run BDSEQ
271 
272  CALL get_code_parms(code_parameters, 'bdseq.xml', '', 'bdseq.xsd')
273  CALL bdseq(eq_in, eq, code_parameters)
274 
275 !... write the results
276 
277 #ifdef MPI
278  CALL mpi_finalize(ierr)
279 #endif
280 
281  stop
282 END PROGRAM bdseq_wrapper
283 
284 
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
bdseq
subroutine bdseq(eq_in, eq_out, code_parameters)
Module implementing a simple equilibrium using CPOs.
Definition: bdseq.F90:6
auxiliaries::prof_init
subroutine prof_init(eq, coreprof)
Definition: wrapper.F90:5
eq
subroutine eq(pcequi, psicon, ncequi, nstep, ngrid,
Definition: Eq_m.f:310
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
bdseq_wrapper
program bdseq_wrapper
Definition: wrapper.F90:208
phys_functions::r_major
real(r8) function r_major(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:102
phys_constants
Definition: phys_constants.f90:1