ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
fusion_dist_sources_module.f90
Go to the documentation of this file.
1 
13 
14 module fusion_dist_sources_module
15 
16  use amns_types
17  integer, parameter :: nreac=4, nspec=6
18  type (amns_handle_type), save :: amns ! AMNS global handle
19  type (amns_handle_rx_type), save :: nuclear(nreac) ! AMNS table handle
20 
21 contains
22 
23  subroutine fusion_dist_sources(coreprof, equilibrium, distsource)
24 
25  use itm_types
26  use amns_module
27  use copy_structures
28  use ets_version
29  use itm_constants
30  use itm_toolbox
31  use distsource_identifier, only: &
32  di_nuclear => nuclear, &
33  di_get_type_name => get_type_name, &
34  di_get_type_description__ind => get_type_description__ind
35 
36  use species_reference_identifier, only: &
37  sr_ion => ion, &
38  sr_electron => electron, &
39  sr_get_type_name => get_type_name, &
40  sr_get_type_description__ind => get_type_description__ind
41 
42  implicit none
43 
44  type (type_coreprof), pointer :: coreprof(:)
45  type (type_equilibrium), pointer :: equilibrium(:)
46  type (type_distsource), pointer :: distsource(:)
47 
48  REAL(R8), DIMENSION(:), POINTER :: rho_eq, rho_core
49 
50  type (amns_reaction_type) :: xx_rx
51  type (amns_query_type) :: query
52  type (amns_answer_type) :: answer
53  type (amns_set_type) :: set
54  type (amns_reactants_type) :: species
55  type (amns_reactant_type) :: nuclei(nspec)
56  real (kind=R8), allocatable :: energy(:), density(:,:), source(:,:), heating(:), rate(:)
57  real (kind=R8), allocatable :: psi(:), volume(:), area(:)
58  integer :: i, j, nr, ns, nucindex, im, iz, nrho, npsi
59  integer, save :: ions_index(nspec), ions_found ! p, D, T, 3He, 4He, n
60  real (kind=R8), save :: mass(nspec) = &
61  (/ itm_mass_h_1, itm_mass_h_2, itm_mass_h_3, itm_mass_he_3, itm_mass_he_4, itm_mn /)
62  integer (kind=R8), save :: rm(4,nreac)
63  character (len=12) :: reaction_label(nreac)
64  real (kind=R8), save :: fusion_energy(nreac)
65  real (kind=R8) :: fraction
66  logical, save :: first=.true.
67  logical :: found, got_psi, got_area, got_volume
68 
69  if(first) then
70  call itm_amns_setup(amns) ! set up the AMNS system
71  query%string = 'version'
72  call itm_amns_query(amns,query,answer)
73  write(*,*) 'fusion: amns data base version = ',trim(answer%string)
74 
75  nuclei(1) = amns_reactant_type(1, 0, 1, 0) ! H
76  nuclei(2) = amns_reactant_type(1, 0, 2, 0) ! D
77  nuclei(3) = amns_reactant_type(1, 0, 3, 0) ! T
78  nuclei(4) = amns_reactant_type(2, 0, 3, 0) ! 3He
79  nuclei(5) = amns_reactant_type(2, 0, 4, 0) ! 4He
80  nuclei(6) = amns_reactant_type(0, 0, 1, 0) ! n
81 
82  rm(:,1) = (/ 2, 2, 1, 3 /) ! D_D_p_T
83  rm(:,2) = (/ 2, 2, 6, 4 /) ! D_D_n_3He
84  rm(:,3) = (/ 2, 3, 6, 5 /) ! D_T_n_4He
85  rm(:,4) = (/ 2, 4, 1, 5 /) ! D_3He_p_4He
86 
87 ! these are the labels defined in "distsource_identifier"
88  reaction_label(1) = 'DD_PT'
89  reaction_label(2) = 'DD_N3He'
90  reaction_label(3) = 'DT_N4He'
91  reaction_label(4) = 'D3He_P4He'
92 
93 ! sanity check for the reactions
94  do nr = 1, nreac
95  im = nint(nuclei(rm(1,nr))%mi) + nint(nuclei(rm(2,nr))%mi) - &
96  nint(nuclei(rm(3,nr))%mi) - nint(nuclei(rm(4,nr))%mi)
97  if(im.ne.0) then
98  write(*,*) 'fusion: baryons not conserved for reaction ', nr
99  stop 'Error'
100  endif
101  iz = nint(nuclei(rm(1,nr))%zn) + nint(nuclei(rm(2,nr))%zn) - &
102  nint(nuclei(rm(3,nr))%zn) - nint(nuclei(rm(4,nr))%zn)
103  if(iz.ne.0) then
104  write(*,*) 'fusion: protons not conserved for reaction ', iz
105  stop 'Error'
106  endif
107  enddo
108 
109 ! calculate the information we will need later for the reactions
110  allocate(species%components(4)) ! set up reactants
111  xx_rx%string='NUC_TT' ! set up reaction
112  xx_rx%isotope_resolved='1'
113 
114  do nr = 1, nreac
115  species%components = (/ &
116  nuclei(rm(1,nr)), &
117  nuclei(rm(2,nr)), &
118  nuclei(rm(3,nr)), &
119  nuclei(rm(4,nr)) /)
120  species%components(3)%LR=1 ! product
121  species%components(4)%LR=1 ! product
122  if(nuclei(rm(3,nr))%zn .eq.0.0_r8) then
123  fraction = mass(rm(3,nr)) / (mass(rm(3,nr))+mass(rm(4,nr)))
124  elseif(nuclei(rm(4,nr))%zn .eq.0.0_r8) then
125  fraction = mass(rm(4,nr)) / (mass(rm(3,nr))+mass(rm(4,nr)))
126  else
127  fraction = 1.0_r8
128  endif
129  fusion_energy(nr) = &
130  (mass(rm(1,nr)) + mass(rm(2,nr)) - mass(rm(3,nr)) - mass(rm(4,nr))) * itm_c**2 * &
131  fraction
132  call itm_amns_setup_table(amns, xx_rx, species, nuclear(nr)) ! set up table
133  query%string = 'reactants'
134  call itm_amns_query_table(nuclear(nr), query, answer)
135  write(*,*) 'fusion: ', nr, trim(answer%string), fusion_energy(nr)/itm_qe/1.0e6
136  enddo
137 
138  deallocate(species%components)
139 
140 ! calculate the mapping vector from external to internal species
141  found = .false.
142  ions_index = 0
143  ions_found = 0
144  do i = 1, size(coreprof(1)%ni%value, dim=2)
145  nucindex = coreprof(1)%compositions%ions(i)%nucindex
146  do ns = 1, nspec
147  found = &
148  (nint(coreprof(1)%compositions%nuclei(nucindex)%zn) .eq. nuclei(ns)%zn) .and. &
149  (nint(coreprof(1)%compositions%nuclei(nucindex)%amn) .eq. nuclei(ns)%mi)
150  if(found) then
151  ions_index(ns) = i
152  ions_found = ions_found+1
153  endif
154  enddo
155  enddo
156  write(*,*) 'fusion: ions_found = ', ions_found
157  write(*,*) 'fusion: ions_index = ', ions_index
158 
159  first = .false.
160 
161  endif
162 
163 ! arrays we will need
164  nrho = size(coreprof(1)%rho_tor)
165  npsi = size(equilibrium(1)%profiles_1d%rho_tor)
166  allocate(energy(nrho))
167  allocate(density(nrho,nspec))
168  allocate(source(nrho,nspec))
169  allocate(heating(nrho))
170  allocate(rate(nrho))
171  allocate(rho_core(nrho))
172  allocate(psi(nrho))
173  allocate(area(nrho))
174  allocate(volume(nrho))
175  allocate(rho_eq(npsi))
176 
177  rho_core = coreprof(1)%rho_tor/coreprof(1)%rho_tor(nrho)
178  rho_eq = equilibrium(1)%profiles_1d%rho_tor/equilibrium(1)%profiles_1d%rho_tor(npsi)
179 
180  got_psi = associated(equilibrium(1)%profiles_1d%psi)
181  if (got_psi) &
182  CALL l3interp(equilibrium(1)%profiles_1d%psi, rho_eq, npsi, &
183  psi, rho_core, nrho)
184  got_area = associated(equilibrium(1)%profiles_1d%area)
185  if (got_area) &
186  CALL l3interp(equilibrium(1)%profiles_1d%area, rho_eq, npsi, &
187  area, rho_core, nrho)
188  got_volume = associated(equilibrium(1)%profiles_1d%volume)
189  if (got_volume) &
190  CALL l3interp(equilibrium(1)%profiles_1d%volume, rho_eq, npsi, &
191  volume, rho_core, nrho)
192 
193 
194 ! energy of the reactants
195  energy = sum(coreprof(1)%ti%value * coreprof(1)%ni%value, 2) / sum(coreprof(1)%ni%value,2)
196 
197 ! densities of reactants
198  density = 0
199  do ns = 1, nspec
200  if(ions_index(ns) .ne. 0) then
201  density(:,ns) = coreprof(1)%ni%value(:,ions_index(ns))
202  endif
203  enddo
204 
205 ! calculate sources
206  source = 0
207  heating = 0
208  do nr = 1, nreac
209  call itm_amns_rx(nuclear(nr), rate, energy) ! get results
210  rate(:) = rate(:) * density(:,rm(1,nr)) * density(:,rm(2,nr))
211  if(rm(1,nr) .eq. rm(2,nr)) rate = rate * 0.5 ! this was missed earlier
212  source(:,rm(1,nr)) = source(:,rm(1,nr)) - rate(:)
213  source(:,rm(2,nr)) = source(:,rm(2,nr)) - rate(:)
214  source(:,rm(3,nr)) = source(:,rm(3,nr)) + rate(:)
215  source(:,rm(4,nr)) = source(:,rm(4,nr)) + rate(:)
216  heating(:) = heating + rate(:) * fusion_energy(nr)
217  enddo
218 
219  write(*,*) 'fusion: density = ', density(1,:)
220  write(*,*) 'fusion: source = ', source(1,:)
221  write(*,*) 'fusion: heating = ', heating(1)
222 
223 ! prepare the output CPO
224  allocate(distsource(1))
225  allocate(distsource(1)%codeparam%codename(1))
226  allocate(distsource(1)%codeparam%codeversion(1))
227  allocate(distsource(1)%codeparam%output_diag(1))
228  distsource(1)%codeparam%codename = 'fusion'
229  distsource(1)%codeparam%codeversion = version
230  distsource(1)%codeparam%output_diag = 'OK'
231  distsource(1)%codeparam%output_flag = 0
232  distsource(1)%time = coreprof(1)%time
233  call copy_cpo(coreprof(1)%compositions, distsource(1)%compositions)
234 
235  allocate(distsource(1)%source(ions_found+1))
236 
237  j = 0
238  do i=1, nspec
239  if(ions_index(i) .ne. 0) then
240  j = j+1
241  allocate(distsource(1)%source(j)%source_id(1))
242  allocate(distsource(1)%source(j)%source_id(1)%type%id(1))
243  allocate(distsource(1)%source(j)%source_id(1)%type%description(1))
244  distsource(1)%source(j)%source_id(1)%type%flag = di_nuclear
245  distsource(1)%source(j)%source_id(1)%type%id = di_get_type_name(di_nuclear)
246  distsource(1)%source(j)%source_id(1)%type%description = di_get_type_description__ind(di_nuclear)
247  allocate(distsource(1)%source(j)%species%type%id(1))
248  allocate(distsource(1)%source(j)%species%type%description(1))
249  distsource(1)%source(j)%species%type%flag = sr_ion
250  distsource(1)%source(j)%species%type%id = sr_get_type_name(sr_ion)
251  distsource(1)%source(j)%species%type%description = sr_get_type_description__ind(sr_ion)
252  distsource(1)%source(j)%species%index = i
253  call copy_cpo(coreprof(1)%rho_tor, distsource(1)%source(j)%profiles_1d%rho_tor)
254  if (got_psi) then
255  allocate(distsource(1)%source(j)%profiles_1d%psi(nrho))
256  distsource(1)%source(j)%profiles_1d%psi = psi
257  endif
258  if (got_area) then
259  allocate(distsource(1)%source(j)%profiles_1d%area(nrho))
260  distsource(1)%source(j)%profiles_1d%area = area
261  endif
262  if (got_volume) then
263  allocate(distsource(1)%source(j)%profiles_1d%volume(nrho))
264  distsource(1)%source(j)%profiles_1d%volume = volume
265  call cubint(nrho, volume, source(:,i), 1, nrho, &
266  distsource(1)%source(j)%global_param%src_rate%value, &
267  distsource(1)%source(j)%global_param%src_rate%abserror)
268  endif
269  allocate(distsource(1)%source(j)%profiles_1d%src_rate%value(nrho))
270  distsource(1)%source(j)%profiles_1d%src_rate%value = source(:,i)
271  endif
272  enddo
273 
274  j = j+1
275  allocate(distsource(1)%source(j)%source_id(1))
276  allocate(distsource(1)%source(j)%source_id(1)%type%id(1))
277  allocate(distsource(1)%source(j)%source_id(1)%type%description(1))
278  distsource(1)%source(j)%source_id(1)%type%flag = di_nuclear
279  distsource(1)%source(j)%source_id(1)%type%id = di_get_type_name(di_nuclear)
280  distsource(1)%source(j)%source_id(1)%type%description = di_get_type_description__ind(di_nuclear)
281  allocate(distsource(1)%source(j)%species%type%id(1))
282  allocate(distsource(1)%source(j)%species%type%description(1))
283  distsource(1)%source(j)%species%type%flag = sr_electron
284  distsource(1)%source(j)%species%type%id = sr_get_type_name(sr_electron)
285  distsource(1)%source(j)%species%type%description = sr_get_type_description__ind(sr_electron)
286  distsource(1)%source(j)%species%index = 0
287  call copy_cpo(coreprof(1)%rho_tor, distsource(1)%source(j)%profiles_1d%rho_tor)
288  if (got_psi) then
289  allocate(distsource(1)%source(j)%profiles_1d%psi(nrho))
290  distsource(1)%source(j)%profiles_1d%psi = psi
291  endif
292  if (got_area) then
293  allocate(distsource(1)%source(j)%profiles_1d%area(nrho))
294  distsource(1)%source(j)%profiles_1d%area = area
295  endif
296  if (got_volume) then
297  allocate(distsource(1)%source(j)%profiles_1d%volume(nrho))
298  distsource(1)%source(j)%profiles_1d%volume = volume
299  call cubint(nrho, volume, heating, 1, nrho, &
300  distsource(1)%source(j)%global_param%src_pow%value, &
301  distsource(1)%source(j)%global_param%src_pow%abserror)
302  endif
303  allocate(distsource(1)%source(j)%profiles_1d%pow_den%value(nrho))
304  distsource(1)%source(j)%profiles_1d%pow_den%value = heating
305 
306 ! deallocate no longer needed internal arrays
307  deallocate(energy, density, source, heating, rate, rho_core, psi, area, volume, rho_eq)
308 
309  return
310 
311  end subroutine fusion_dist_sources
312 
313  subroutine fusion_finalize()
314 
315  use amns_module
316 
317  implicit none
318 
319  integer :: nr
320 
321  do nr = 1, nreac
322  call itm_amns_finish_table(nuclear(nr)) ! finish with table
323  enddo
324  call itm_amns_finish(amns) ! finish with amns
325 
326  end subroutine fusion_finalize
327 
328 end module fusion_dist_sources_module
fusion_dist_sources_module
Module implementing fusion sources.
Definition: fusion_dist_sources_module.f90:14
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
fusion_dist_sources_module::fusion_dist_sources
subroutine fusion_dist_sources(coreprof, equilibrium, distsource)
Definition: fusion_dist_sources_module.f90:23
itm_toolbox
prototype ITM_TOOLBOX
Definition: itm_toolbox.F90:11
fusion_dist_sources_module::fusion_finalize
subroutine fusion_finalize()
Definition: fusion_dist_sources_module.f90:313
ets_version
Definition: ets_version.f90:1
cubint
subroutine cubint(ntab, xtab, ftab, ia_in, ib_in, result, error)
Definition: cubint.f90:1