Commit eeadae21 authored by Fabian Wieschollek's avatar Fabian Wieschollek

First files

parents
!> particles/examples/W_example.f90
!> Requires: acd50_w.dat, scd50_w.dat and the JOREK input file
!> Compile with `j2p W_example` and run with `W_example < in1`
!> Read fields from jorek00721.h5 and onwards and push W particles
program Radiative
use mod_openadas
use constants
use mumps_module
use pastix_module
use wsmp_module
use data_structure
use phys_module
use mod_parameters
use mod_log_params
use global_distributed_matrix
use nodes_elements
use pellet_module
use equil_info
use mod_boundary, only: boundary_from_grid
use vacuum
use vacuum_response, only: get_vacuum_response, update_response, init_wall_currents, I_coils
use vacuum_equilibrium, only: import_external_fields
use live_data
use mod_bootstrap_functions
use construct_matrix_mod, only : construct_matrix
use mod_global_matrix_structure
use mod_import_restart
use mod_export_restart
use mod_element_rtree, only: populate_element_rtree
use mod_interp
use basis_at_gaussian, only: initialise_basis
use mod_expression, only: exprs_all_int, init_expr
use mod_integrals3D
use mod_interp_splinear
use mod_non_coronal
use mod_coronal
use mpi_mod
implicit none
type(adf11_all) :: adas
type(non_coronal) :: ncor
character(len=128) :: filename
character(len=23), dimension(3) :: header = (/ 't','T', 'P_rad' /)
real*8 :: dtime, stime, ftime, tau_diff
integer :: step, nlines
! Conversion: Time JOREK->SI : 1.d0 = sqrt_m0_rho0 s
! Temperature JOREK->SI : 1.d0 = 1/(2.d0*K_BOLTZ*MU_ZERO*central_density*1.d20) K
! SI->eV : 1.d0 = 8.61733326214517743366d-5 eV
! Density JOREK->SI : 1.d0 = 1.d20*central_density m^-3
! central_density = 1.d0
! sqrt_m0_rho0 = 6.483638667465844d-007
! 1/(2.d0*K_BOLTZ*MU_ZERO*1.d20) = 288188681.82099485
! --- Prepare the evolution
write(*,*) '#### Prepare Evolution ####'
adas = read_adf11('89_w') ! ADAS Data for choosen Impurity
ncor = non_coronal(adas,20.d0,4.064627111056752d0) ! Non-Coronal-Type with initial n_e=1.d0 and T=1.e-3
stime = 0.d0 ! Simulation time (s)
tau_diff = 3.d99 ! Diffusion time (s)
filename = "tungstentest-SI.dat" ! Name of outputfile
nlines = 1 ! output every ...
dtime = 1.d0 ! 0.01d0 ! Timestep
ftime = 10000.d0 !10000.d0 ! Final time
! --- Output first line(s) for stimet=0
open(20,file=filename)
write(20,'(A,3A23)') '# ',header
call ncor%output(20)
! --- Start the Evolution
write(*,*) '##### Start Evolution #####'
write(*,*) '##### #Steps:',int(ftime/dtime),'####'
write(*,*) '##### #Lines:',int(ftime/dtime)/nlines,'###'
do step=1,int(ftime/dtime)
! --- Time step
stime=step*dtime
call ncor%update(dtime,tau_diff)
! --- Output
if ( int(mod(step,nlines)) .eq. 0 ) call ncor%output(20)
! --- Update T
call ncor%update_T(log10(10.d0**4.064627111056752d0*(1.d0+step)))
call ncor%fractions_coronal()
end do
close (20)
end program Radiative
!> particles/examples/W_example.f90
!> Requires: acd50_w.dat, scd50_w.dat and the JOREK input file
!> Compile with `j2p W_example` and run with `W_example < in1`
!> Read fields from jorek00721.h5 and onwards and push W particles
program Radiative
use mod_openadas
use constants
use mumps_module
use pastix_module
use wsmp_module
use data_structure
use phys_module
use mod_parameters
use mod_log_params
use global_distributed_matrix
use nodes_elements
use pellet_module
use equil_info
use mod_boundary, only: boundary_from_grid
use vacuum
use vacuum_response, only: get_vacuum_response, update_response, init_wall_currents, I_coils
use vacuum_equilibrium, only: import_external_fields
use live_data
use mod_bootstrap_functions
use construct_matrix_mod, only : construct_matrix
use mod_global_matrix_structure
use mod_import_restart
use mod_export_restart
use mod_element_rtree, only: populate_element_rtree
use mod_interp
use basis_at_gaussian, only: initialise_basis
use mod_expression, only: exprs_all_int, init_expr
use mod_integrals3D
use mod_interp_splinear
use mod_non_coronal
use mod_coronal
use mpi_mod
implicit none
type(adf11_all) :: adas
type(non_coronal) :: ncor
character(len=128) :: filename
character(len=23), dimension(7) :: header = (/ 't','log10_Te','log10_ne','Li0','Li1','Li2','Li3' /)
real*8 :: dtime, stime, ftime, tau_diff
integer :: step, nlines
! Conversion: Time JOREK->SI : 1.d0 = sqrt_m0_rho0 s
! Temperature JOREK->SI : 1.d0 = 1/(2.d0*K_BOLTZ*MU_ZERO*central_density*1.d20) K
! Density JOREK->SI : 1.d0 = 1.d20*central_density m^-3
! central_density = 1.d0
! sqrt_m0_rho0 = 6.483638667465844d-007
! 1/(2.d0*K_BOLTZ*MU_ZERO*1.d20) = 288188681.82099485
! --- Prepare the evolution
write(*,*) '#### Prepare Evolution ####'
adas = read_adf11('96_li') ! ADAS Data for choosen Impurity
ncor = non_coronal(adas,20.d0,5.459676920547336d0) ! Non-Coronal-Type with initial n_e=1.d0 and T=1.e-3
stime = 0.d0 ! Simulation time (s)
tau_diff = 3.d99 ! Diffusion time (s)
filename = "lithiumtest-SI.dat" ! Name of outputfile
nlines = 100 ! output every ...
dtime = 6.483638667465844d-09 ! 0.01d0 ! Timestep
ftime = 0.006483638667465844d0 !10000.d0 ! Final time
! --- Output first line(s) for stimet=0
open(20,file=filename)
write(20,'(A,13A23)') '# ',header
call ncor%output(20)
! --- Start the Evolution
write(*,*) '##### Start Evolution #####'
write(*,*) '##### #Steps:',int(ftime/dtime),'####'
write(*,*) '##### #Lines:',int(ftime/dtime)/nlines,'###'
do step=1,int(ftime/dtime)
! --- Time step
stime=step*dtime
call ncor%update(dtime,tau_diff)
! --- Output
if ( int(mod(step,nlines)) .eq. 0 ) call ncor%output(20)
! --- Update T
!call ncor%update_T(1.d-6*(1.d0+step))
end do
close (20)
end program Radiative
!> module takes the OPEN-ADAS data to calculate the non_coronal equil-
!> ibrium temperature and radiation.
!> If you need time-dependent solutions of the corona matrix timestepping look
!> in the revision history for this file.
module mod_non_coronal
use mod_openadas
use mod_interp_splinear
use constants
implicit none
private
public non_coronal
! -------------------------------------------------------------------------------------
! --- Definition of the datatype
! -------------------------------------------------------------------------------------
type non_coronal
type (ADF11_all) :: ad !< ADF11 datatype
real*8, allocatable :: fractions(:) !< Fractional charge states. Should sum to 1 but we do not check it!
real*8, allocatable :: coronal(:) !< charge states regarding coronal equilibrium. Should sum to 1 but we do not check it!
real*8, allocatable :: ion_rate(:) !< Current ionization rates (1/s)
real*8, allocatable :: rec_rate(:) !< Current recombinations rates (1/s)
real*8, allocatable :: rad(:) !< Current radiation coefficients
real*8, allocatable :: Z(:) !< The charge number at each charge state
real*8 :: density !< current log10 density (m^-3)
real*8 :: temperature !< current log10 temperature (K)
real*8 :: time !< current time (s)
real*8 :: Prad !< current log10 Radiated power per ion (W)
integer :: n_Z0 !< Z plus 1 to factor in the atomic state
contains
procedure :: update => update_fractions
procedure :: fractions_coronal => fractions_coronal
procedure :: update_T => update_T
procedure :: update_N => update_N
procedure :: update_GCRs => update_GCRs
procedure :: update_c => update_coronal
procedure :: output => output_non_coronal
procedure :: non_coronal_Prad => non_coronal_Prad
procedure :: z_EFF => z_EFF
end type non_coronal
interface non_coronal
module procedure non_coronal_equilibrium
end interface non_coronal
contains
! -------------------------------------------------------------------------------------
! --- Setup the datatype
! -------------------------------------------------------------------------------------
function non_coronal_equilibrium(ad,density,temperature) result(cor)
use constants
type (ADF11_all), intent(in) :: ad !< ADF11 datatype
real*8, intent(in) :: density !< initial log10 density in m^-3
real*8, intent(in) :: temperature !< initial log10 electron temperature in K
type (non_coronal) :: cor !< non_coronal equilibrium datatype
integer :: iz
cor%ad = ad
cor%n_Z0 = ad%n_Z+1
allocate(cor%fractions(cor%n_Z0),cor%coronal(cor%n_Z0),cor%Z(cor%n_Z0),cor%ion_rate(ad%n_Z),cor%rec_rate(ad%n_Z),cor%rad(ad%n_Z))
do iz=1,cor%n_Z0
cor%Z(iz) = real(iz,8)
enddo
cor%time = 0.d0
cor%fractions = 0.0d0
cor%fractions(1) = 1.0d0
cor%density = density
cor%temperature = temperature
call cor%update_GCRs()
end function non_coronal_equilibrium
! -------------------------------------------------------------------------------------
! --- Update Temperature
! -------------------------------------------------------------------------------------
subroutine update_T(cor,temperature)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: temperature !< log10 temperature in m^-3
cor%temperature=temperature
call cor%update_GCRs()
end subroutine update_T
! -------------------------------------------------------------------------------------
! --- Update Density
! -------------------------------------------------------------------------------------
subroutine update_N(cor,density)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: density !< log10 density in m^-3
cor%density=density
call cor%update_GCRs()
end subroutine update_N
! -------------------------------------------------------------------------------------
! --- Update Coefficients
! -------------------------------------------------------------------------------------
subroutine update_GCRs(cor)
class(non_coronal), intent(inout) :: cor
integer :: iz
real*8 :: rad_RB, rad_LT
do iz=1,cor%ad%n_Z
call cor%ad%ACD%interp(iz, cor%density, cor%temperature, cor%rec_rate(iz))
call cor%ad%SCD%interp(iz, cor%density, cor%temperature, cor%ion_rate(iz))
call cor%ad%PRB%interp(iz+1, cor%density, cor%temperature, rad_RB)
call cor%ad%PLT%interp(iz+1, cor%density, cor%temperature, rad_LT)
cor%rad(iz) = rad_RB + rad_LT
end do
! From m^3/s to 1/s
cor%rec_rate=cor%rec_rate*10**cor%density
cor%ion_rate=cor%ion_rate*10**cor%density
call cor%update_c()
end subroutine update_GCRs
! -------------------------------------------------------------------------------------
! --- Update f coronal
! -------------------------------------------------------------------------------------
subroutine update_coronal(cor)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
integer :: iz
cor%coronal = 0.d0
cor%coronal(1) = 1.0d0
do iz=2,cor%n_Z0
cor%coronal(iz) = cor%coronal(iz-1) * cor%ion_rate(iz-1)/cor%rec_rate(iz-1)
end do
cor%coronal = cor%coronal / sum(cor%coronal)
end subroutine update_coronal
! -------------------------------------------------------------------------------------
! --- Update f temporal evolution
! -------------------------------------------------------------------------------------
subroutine update_fractions(cor,dtime,tau)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: dtime, tau !< time increment in s, addition impurites
real*8, dimension(cor%n_Z0) :: fractionsN
integer :: iz
cor%time = cor%time+dtime
fractionsN = cor%fractions
do iz=1,cor%n_Z0
if ( iz .gt. 1) then
fractionsN(iz) = fractionsN(iz) &
+ cor%fractions(iz-1) * cor%ion_rate(iz-1) * dtime & ! Gain by Ionization from (iz-1) to iz
- cor%fractions(iz) * cor%rec_rate(iz-1) * dtime ! Loss by Recombination from iz to (iz-1)
end if
if ( iz .lt. cor%n_Z0) then
fractionsN(iz) = fractionsN(iz) &
+ cor%fractions(iz+1) * cor%rec_rate(iz) * dtime & ! Gain by Recombination from (iz+1) to iz
- cor%fractions(iz) * cor%ion_rate(iz) * dtime ! Loss by Ionization from iz to (iz+1)
end if
end do
cor%fractions = fractionsN
! Diffusion and fueling
!cor%fractions = cor%fractions*(1.-dtime/tau)
!cor%fractions(1) = (1.-sum(cor%fractions(2:cor%n_Z0)))
end subroutine update_fractions
! -------------------------------------------------------------------------------------
! --- Set f temporal to f coronal
! -------------------------------------------------------------------------------------
subroutine fractions_coronal(cor)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
integer :: iz
call cor%update_c()
do iz=1,cor%n_Z0
cor%fractions(iz) = cor%coronal(iz)
end do
end subroutine fractions_coronal
! -------------------------------------------------------------------------------------
! --- Calculates P_rad for given f
! -------------------------------------------------------------------------------------
function non_coronal_Prad(cor,fractions)
class(non_coronal), intent(in) :: cor
real*8, dimension(cor%n_Z0) :: fractions
real*8 :: non_coronal_Prad !< Output power in W / atom
non_coronal_Prad = dot_product(fractions(2:cor%n_Z0), cor%rad)
end function non_coronal_Prad
! -------------------------------------------------------------------------------------
! --- Calculates z_EFF for given f
! -------------------------------------------------------------------------------------
function z_EFF(cor, fractions)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
real*8, dimension(cor%n_Z0) :: fractions !< The charge number at each charge state
real*8 :: z_EFF
z_EFF = dot_product(fractions,cor%Z)
end function z_EFF
! ----------------------------------------------------------------------------------
! --- Print data out
! ----------------------------------------------------------------------------------
subroutine output_non_coronal(cor,pointer)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
integer, intent(in) :: pointer
integer :: iz
real*8 :: Prad_NC, Prad_C , zeff_NC, zeff_C
! --- 0 stime
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%time
! --- 1 temperature
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%temperature
! --- 2 Prad
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%non_coronal_Prad(cor%coronal)
! --- new line
write(pointer,'(a)') ' '
end subroutine output_non_coronal
end module mod_non_coronal
This diff is collapsed.
!> module takes the OPEN-ADAS data to calculate the non_coronal equil-
!> ibrium temperature and radiation.
!> If you need time-dependent solutions of the corona matrix timestepping look
!> in the revision history for this file.
module mod_non_coronal_si
use mod_openadas
use mod_interp_splinear
use constants
implicit none
private
public non_coronal
! -------------------------------------------------------------------------------------
! --- Definition of the datatype
! -------------------------------------------------------------------------------------
type non_coronal
type (ADF11_all) :: ad !< ADF11 datatype
real*8, allocatable :: fractions(:) !< Fractional charge states. Should sum to 1 but we do not check it!
real*8, allocatable :: coronal(:) !< charge states regarding coronal equilibrium. Should sum to 1 but we do not check it!
real*8, allocatable :: ion_rate(:) !< Current ionization rates (1/s)
real*8, allocatable :: rec_rate(:) !< Current recombinations rates (1/s)
real*8, allocatable :: rad(:) !< Current radiation coefficients
real*8, allocatable :: Z(:) !< The charge number at each charge state
real*8 :: density !< current log10 density (m^-3)
real*8 :: temperature !< current log10 temperature (K)
real*8 :: time !< current time (s)
real*8 :: Prad !< current log10 Radiated power per ion (W)
integer :: n_Z0 !< Z plus 1 to factor in the atomic state
contains
procedure :: update => update_fractions
procedure :: fractions_coronal => fractions_coronal
procedure :: update_T => update_T
procedure :: update_N => update_N
procedure :: update_GCRs => update_GCRs
procedure :: update_c => update_coronal
procedure :: output => output_non_coronal
procedure :: non_coronal_Prad => non_coronal_Prad
procedure :: z_EFF => z_EFF
end type non_coronal
interface non_coronal
module procedure non_coronal_equilibrium
end interface non_coronal
contains
! -------------------------------------------------------------------------------------
! --- Setup the datatype
! -------------------------------------------------------------------------------------
function non_coronal_equilibrium(ad,density,temperature) result(cor)
use constants
type (ADF11_all), intent(in) :: ad !< ADF11 datatype
real*8, intent(in) :: density !< initial log10 density in m^-3
real*8, intent(in) :: temperature !< initial log10 electron temperature in K
type (non_coronal) :: cor !< non_coronal equilibrium datatype
integer :: iz
cor%ad = ad
cor%n_Z0 = ad%n_Z+1
allocate(cor%fractions(cor%n_Z0),cor%coronal(cor%n_Z0),cor%Z(cor%n_Z0),cor%ion_rate(ad%n_Z),cor%rec_rate(ad%n_Z),cor%rad(ad%n_Z))
do iz=1,cor%n_Z0
cor%Z(iz) = real(iz,8)
enddo
cor%time = 0.d0
cor%fractions = 0.0d0
cor%fractions(1) = 1.0d0
cor%density = density
cor%temperature = temperature
call cor%update_GCRs()
end function non_coronal_equilibrium
! -------------------------------------------------------------------------------------
! --- Update Temperature
! -------------------------------------------------------------------------------------
subroutine update_T(cor,temperature)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: temperature !< log10 temperature in m^-3
cor%temperature=temperature
call cor%update_GCRs()
end subroutine update_T
! -------------------------------------------------------------------------------------
! --- Update Density
! -------------------------------------------------------------------------------------
subroutine update_N(cor,density)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: density !< log10 density in m^-3
cor%density=density
call cor%update_GCRs()
end subroutine update_N
! -------------------------------------------------------------------------------------
! --- Update Coefficients
! -------------------------------------------------------------------------------------
subroutine update_GCRs(cor)
class(non_coronal), intent(inout) :: cor
integer :: iz
real*8 :: rad_RB, rad_LT
do iz=1,cor%ad%n_Z
call cor%ad%ACD%interp(iz, cor%density, cor%temperature, cor%rec_rate(iz))
call cor%ad%SCD%interp(iz, cor%density, cor%temperature, cor%ion_rate(iz))
call cor%ad%PRB%interp(iz+1, cor%density, cor%temperature, rad_RB)
call cor%ad%PLT%interp(iz+1, cor%density, cor%temperature, rad_LT)
cor%rad(iz) = rad_RB + rad_LT
end do
! From m^3/s to 1/s
cor%rec_rate=cor%rec_rate*10**cor%density
cor%ion_rate=cor%ion_rate*10**cor%density
call cor%update_c()
end subroutine update_GCRs
! -------------------------------------------------------------------------------------
! --- Update f coronal
! -------------------------------------------------------------------------------------
subroutine update_coronal(cor)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
integer :: iz
cor%coronal = 0.d0
cor%coronal(1) = 1.0d0
do iz=2,cor%n_Z0
cor%coronal(iz) = cor%coronal(iz-1) * cor%ion_rate(iz-1)/cor%rec_rate(iz-1)
end do
cor%coronal = cor%coronal / sum(cor%coronal)
end subroutine update_coronal
! -------------------------------------------------------------------------------------
! --- Update f temporal evolution
! -------------------------------------------------------------------------------------
subroutine update_fractions(cor,dtime,tau)
class(non_coronal), intent(inout) :: cor !< Coronal equilibrium type
real*8, intent(in) :: dtime, tau !< time increment in s, addition impurites
real*8, dimension(cor%n_Z0) :: fractionsN
integer :: iz
cor%time = cor%time+dtime
fractionsN = cor%fractions
do iz=1,cor%n_Z0
if ( iz .gt. 1) then
fractionsN(iz) = fractionsN(iz) &
+ cor%fractions(iz-1) * cor%ion_rate(iz-1) * dtime & ! Gain by Ionization from (iz-1) to iz
- cor%fractions(iz) * cor%rec_rate(iz-1) * dtime ! Loss by Recombination from iz to (iz-1)
end if
if ( iz .lt. cor%n_Z0) then
fractionsN(iz) = fractionsN(iz) &
+ cor%fractions(iz+1) * cor%rec_rate(iz) * dtime & ! Gain by Recombination from (iz+1) to iz
- cor%fractions(iz) * cor%ion_rate(iz) * dtime ! Loss by Ionization from iz to (iz+1)
end if
end do
cor%fractions = fractionsN
! Diffusion and fueling
!cor%fractions = cor%fractions*(1.-dtime/tau)
!cor%fractions(1) = (1.-sum(cor%fractions(2:cor%n_Z0)))
end subroutine update_fractions
! -------------------------------------------------------------------------------------
! --- Set f temporal to f coronal
! -------------------------------------------------------------------------------------