Commit d76b5e04 authored by Fabian Wieschollek's avatar Fabian Wieschollek

Greater Changes!

parent 6ebc5841
&input
ADASname = '89_ar' ! Argon
outputname = 'CE-argon.dat'
output_commands = 'Te', 'fractions', 'L_rad', '<Z>', 'Z_eff' !Temperature for x-axis
line_zero = .f.
unit = 1
nout = 1
tstep = 1.0d-1 ! Go in timesteps of 0.1 i.e. 0.1 eV
ftime = 10000 ! Final time, where Te=10keV
temperature = 2.013320969d-5 ! Begin with 1 eV initially
density_plasma = 1.d0 ! For ne=1.d20
density_imp = 1.d1
t_rad = 1.d99 ! Turn off radiative cooling terms
t_equi = -1.d99 ! Assume coronal equilibrium from the beginning
t_pres = 1.d99 ! Turn off adiabatic terms, keep densities
C_cooling(1) = -2.013320969d-5 ! Increase temperature by 1 eV per time unit
C_cooling(2) = 0.d0 ! Start cooling from the beginning
C_cooling(3) = 1.d-5 ! Stop cooling, if Te < ...
/
#####################################################
+++++++ OPTIONS FOR OUTPUT_COMMANDS (UP TO 10) ++++++
#####################################################
't' > time
'Te' > temperature
'ne' > electron density
'fractions' > fractions of imp*
'ion_rate' > ionization rates*
'ion_rate' > recomb rates*
'Z_eff' > effectice Z
'<Z>' > avergaed Z
'L_rad' > radiation rate
'P_rad' > radiated P density
'log10_Te' > log10(Te/K)
'log10_ne' > log10(ne/m3)
* creates Z+1 columns for each charge state!
#####################################################
&input
outputname = 'lithiumtest.dat' ! Name of outputfile
ADASname = '96_li' ! Prefix of ADAS files
ADASname = '96_li' ! Prefix of ADAS files
tstep = 0.01d0 ! Timestep
ftime = 10000.d0 ! Final time
nout = 100 ! Output every
outputname = 'Li-Example.dat' ! Name of outputfile
output_commands = 't', 'Te', 'ne' 'fractions' ! output commands, see commands.README
line_zero = .t. ! Print line for t=0 with initial conidtions?
unit = 1 ! 0: JOREK-units 1: SI-units with eV for Th 2: SI-units with K for Th
nout = 100 ! Output every
temperature = 0.04d0 ! Initial electron temperature (everything in JOREK-units, central_density=1.d0)
density_plasma = 1.d0 ! Plasma ion density
density_imp = 1.d-1 ! Impurity ion density
tstep = 0.01d0 ! Timestep
ftime = 10000.d0 ! Final time
t_rad = 0.d00 ! Turn on radiative cooling term after
t_equi = 1.d99 ! Assume coronal equilibrium after
temperature = 0.04d0 ! Initial electron temperature (everything in JOREK-units, central_density=1.d0)
density_plasma = 1.d0 ! Plasma ion density
density_imp = 1.d-1 ! Impurity ion density
C_cooling(1) = 9.d-6 ! Linear Cooling rate dT/dt
C_cooling(2) = 1.d3 ! Start cooling after ...
C_cooling(3) = 1.d-2 ! Stop cooling, if Te < ...
t_rad = -1.d00 ! Turn on radiative cooling term after
t_equi = 1.d99 ! Assume coronal equilibrium after
t_pres = -1.d00 ! Assume adiabatic, evolve densities consistently after
C_cooling(1) = 9.d-6 ! Linear Cooling rate dT/dt
C_cooling(2) = 1.d3 ! Start cooling after ...
C_cooling(3) = 1.d-2 ! Stop cooling, if Te < ...
/
......@@ -8,28 +8,31 @@ program Radiative
type(non_coronal) :: ncor
character(len=128) :: outputname
character(len=5) :: ADASname
real*8 :: tstep, ftime, temperature, density_plasma, density_imp, t_rad, t_equi
real*8 :: tstep, ftime, temperature, density_plasma, density_imp, t_rad, t_equi, t_pres
real*8, dimension(3) :: C_cooling
integer :: step, nout, ierr
integer :: step, nout, ierr, unit
logical :: line_zero
character(len=16), dimension(OUTPUT_MAX) :: output_commands
namelist /input/ outputname, ADASname, &
namelist /input/ outputname, output_commands, ADASname, &
tstep, ftime, temperature, density_plasma, density_imp, t_rad, t_equi, &
C_cooling, nout
C_cooling, nout, line_zero, t_pres, unit
! --- Prepare the evolution
write(*,*) '#### Prepare Evolution #####'
read(5,input)
adas = read_adf11(ADASname) ! ADAS Data for choosen Impurity
ncor = non_coronal(adas,temperature,density_plasma,density_imp,t_rad,t_equi,C_cooling) ! Non-Coronal-Type
adas = read_adf11(ADASname) ! ADAS Data for choosen Impurity
ncor = non_coronal(adas,temperature,density_plasma,density_imp,t_rad,t_equi,t_pres,C_cooling) ! Non-Coronal-Type
! --- Output first line for stime = 0
open(20,file=outputname)
call ncor%output_header(20)
call ncor%output(20)
call ncor%output_header(20,output_commands)
if ( line_zero .eq. .true. ) then
call ncor%output(20,output_commands,unit)
end if
......@@ -38,11 +41,12 @@ write(*,*) '##### Start Evolution #####'
write(*,*) '##### #Steps:',int(ftime/tstep),'####'
write(*,*) '##### #Lines:',int(ftime/tstep)/nout,'###'
do step=1,int(ftime/tstep)
! --- Time step
call ncor%update(tstep)
! --- Output
if ( int(mod(step,nout)) .eq. 0 ) call ncor%output(20)
if ( int(mod(step,nout)) .eq. 0 ) call ncor%output(20,output_commands,unit)
end do
......
......@@ -7,7 +7,7 @@ use mod_openadas
use mod_interp_splinear
use constants
implicit none
private
public non_coronal
! -------------------------------------------------------------------------------------
......@@ -17,23 +17,26 @@ type non_coronal
type (ADF11_all) :: ad !< ADF11 datatype
!Simulation variables
real*8 :: time !< current time
real*8 :: temperature !< current electron temperature
real*8 :: density !< current electrion density
real*8 :: density_plasma !< constant density of ions of the background plasma
real*8 :: density_imp !< constant density of the imuprity ions
real*8, dimension(3) :: C_cooling !< Coefficient for linear aritifical cooling C(1): dT/dt C(2): when starts cooling, C(3): Tmin
real*8 :: t_rad, t_equi !< when starts radiative cooling, when go into coronal equilbrium
real*8, allocatable :: fractions(:) !< Fractional charge states of impurites
real*8 :: time !< current time
real*8 :: temperature !< current electron temperature
real*8 :: density !< current electrion density
real*8 :: density_plasma !< constant density of ions of the background plasma
real*8 :: density_imp !< constant density of the imuprity ions
real*8, dimension(3) :: C_cooling !< Coefficient for linear aritifical cooling C(1): dT/dt C(2): when starts cooling, C(3): Tmin
real*8 :: t_rad, t_equi, t_pres !< when starts radiative cooling, when go into coronal equilbrium, when starts adiabatic
real*8, allocatable :: fractions(:) !< Fractional charge states of impurites
! Constants
integer :: n_Z0 !< Z plus 1 to factor in the atomic state
real*8, allocatable :: Z(:) !< The charge number at each charge state
real*8, allocatable :: ZZ(:) !< The squared charge number at each charge state
! Temporary auxiliary values
real*8, allocatable :: ion_rate(:) !< Current ionization rates
real*8, allocatable :: rec_rate(:) !< Current recombinations rates
real*8, allocatable :: rad(:) !< Current radiation coefficients
real*8, allocatable :: rad_LT(:) !< Current radiation coefficients
real*8, allocatable :: rad_RB(:) !< Current radiation coefficients
contains
procedure :: set_T => set_T
......@@ -45,7 +48,9 @@ contains
procedure :: calculate_GCRs => calculate_GCRs
procedure :: Prad => calculate_Prad
procedure :: Lrad => calculate_Lrad
procedure :: z_EFF => calculate_z_EFF
procedure :: z_avg => calculate_z_avg
procedure :: equilibrium => calculate_equilibrium
procedure :: cooling => calculate_cooling
end type non_coronal
......@@ -59,25 +64,26 @@ contains
! -------------------------------------------------------------------------------------
! --- Setup the datatype
! -------------------------------------------------------------------------------------
function non_coronal_equilibrium(ad,temperature,density_plasma,density_imp,t_rad,t_equi,C_cooling) result(cor)
function non_coronal_equilibrium(ad,temperature,density_plasma,density_imp,t_rad,t_equi,t_pres,C_cooling) result(cor)
use constants
type (ADF11_all), intent(in) :: ad !< ADF11 datatype
real*8, intent(in) :: density_plasma !< plasma (ion) density
real*8, intent(in) :: density_imp !< impurity (ion) density
real*8, intent(in) :: temperature !< initial temperature
real*8, intent(in) :: t_rad, t_equi !< when starts radiative cooling, when go into coronal equilbrium
real*8, dimension(3), intent(in) :: C_cooling !< Coefficient for linear aritifical cooling C(1) dT/dt C(2) when starts cooling, C(3) Tmin
!real*8, intent(in) :: tau_diff !< Characteristic Diffusion time
type (non_coronal) :: cor !< non_coronal equilibrium datatype
type (ADF11_all), intent(in) :: ad !< ADF11 datatype
real*8, intent(in) :: density_plasma !< plasma (ion) density
real*8, intent(in) :: density_imp !< impurity (ion) density
real*8, intent(in) :: temperature !< initial temperature
real*8, intent(in) :: t_rad, t_equi, t_pres !< when starts radiative cooling, when go into coronal equilbrium
real*8, dimension(3), intent(in) :: C_cooling !< Coefficient for linear aritifical cooling C(1) dT/dt C(2) when starts cooling, C(3) Tmin
!real*8, intent(in) :: tau_diff !< Characteristic Diffusion time
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%Z(cor%n_Z0),cor%ion_rate(ad%n_Z),cor%rec_rate(ad%n_Z),cor%rad(ad%n_Z))
allocate(cor%fractions(cor%n_Z0),cor%Z(cor%n_Z0),cor%ZZ(cor%n_Z0),cor%ion_rate(ad%n_Z),cor%rec_rate(ad%n_Z),cor%rad(ad%n_Z),cor%rad_LT(ad%n_Z),cor%rad_RB(ad%n_Z))
do iz=1,cor%n_Z0
cor%Z(iz) = real(iz-1,8)
cor%ZZ(iz) = real(iz-1,8)**2
enddo
cor%time = 0.d0
......@@ -90,6 +96,7 @@ cor%fractions(1) = 1.0d0
cor%C_cooling = C_cooling
cor%t_rad = t_rad
cor%t_equi = t_equi
cor%t_pres = t_pres
call cor%calculate_GCRs()
......@@ -200,6 +207,10 @@ else
ddensity_dt = cor%density_imp * ( dot_product(cor%fractions,cor%Z) - dot_product(fractions0,cor%Z) ) / dtime
end if
if ( cor%time .lt. cor%t_pres ) then
ddensity_dt = 0.d0
end if
!Diffusion and fueling
!dfractions_dt = cor%fractions/tau
!dfractions_dt(1) = sum(cor%dfractions_dt(2:cor%n_Z0))
......@@ -215,6 +226,12 @@ cor%density = density0 + ddensity_dt * dtime
cor%time = cor%time + dtime
cor%temperature = temperature0 + dtemperature_dt * dtime
do iz=1,cor%n_Z0
if ( cor%fractions(iz) .lt. 1.d-99 ) then
cor%fractions(iz) = 0.d0
end if
end do
call cor%calculate_GCRs()
end subroutine update_main
......@@ -222,27 +239,51 @@ end subroutine update_main
! ----------------------------------------------------------------------------------
! --- Print header
! ----------------------------------------------------------------------------------
subroutine output_header(cor,pointer)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
integer, intent(in) :: pointer
integer :: iz
character(len=21) :: header_Z
subroutine output_header(cor,pointer,output_commands)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
integer, intent(in) :: pointer
character(len=16), dimension(OUTPUT_MAX), intent(in) :: output_commands
integer :: iz, c
character(len=21) :: header_Z
write(pointer,'(A)',advance='no') '# '
! --- 0 stime
write(pointer,'(A23)',advance='no') 'time'
! --- 1 temperature
write(pointer,'(A23)',advance='no') 'temperature'
! --- 2 temperature
write(pointer,'(A23)',advance='no') 'density'
! --- 2+n_Z0 to 3+2*n_Z0-1 f states
do iz = 1, cor%n_Z0
write(header_Z,'(A10,I2,A1)') 'fractions(',iz,')'
write(pointer,'(A23)',advance='no') header_Z
do c=1,OUTPUT_MAX
select case( output_commands(c) )
case ('t')
write(pointer,'(A21,A2)',advance='no') 'time',' '
case ('Te')
write(pointer,'(A21,A2)',advance='no') 'temperature',' '
case ('ne')
write(pointer,'(A23)',advance='no') 'density '
case ('fractions')
do iz = 1, cor%n_Z0
write(header_Z,'(A10,I2,A1)') 'fractions(',iz,') '
write(pointer,'(A23)',advance='no') header_Z
end do
case ('ion_rate')
do iz = 1, cor%ad%n_Z
write(header_Z,'(A10,I2,A1)') 'ion_rate(',iz,') '
write(pointer,'(A23)',advance='no') header_Z
end do
case ('rec_rate')
do iz = 1, cor%ad%n_Z
write(header_Z,'(A10,I2,A1)') 'rec_rate(',iz,') '
write(pointer,'(A23)',advance='no') header_Z
end do
case ('Z_eff')
write(pointer,'(A23)',advance='no') 'Z_eff '
case ('<Z>')
write(pointer,'(A23)',advance='no') '<Z> '
case ('P_rad')
write(pointer,'(A23)',advance='no') 'P_rad '
case ('L_rad')
write(pointer,'(A23)',advance='no') 'L_rad '
case ('log10_Te')
write(pointer,'(A23)',advance='no') 'log10_Te '
case ('log10_ne')
write(pointer,'(A23)',advance='no') 'log10_ne '
end select
end do
! --- new line
......@@ -254,26 +295,77 @@ end subroutine output_header
! ----------------------------------------------------------------------------------
! --- Print data out for a time step
! ----------------------------------------------------------------------------------
subroutine output_main(cor,pointer)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
integer, intent(in) :: pointer
integer :: iz
subroutine output_main(cor,pointer,output_commands,unit)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
integer, intent(in) :: pointer
character(len=16), dimension(OUTPUT_MAX), intent(in) :: output_commands
integer, intent(in) :: unit
integer :: iz, c
real*8 :: z_avg, z_EFF, P_rad, log10_Te, log10_ne, L_rad
real*8 :: norm_Th, norm_T, norm_M, norm_N, norm_E
if ( unit .lt. 1 ) then
norm_Th = 1.d0
norm_T = 1.d0
norm_N = 1.d0
norm_M = 1.d0
norm_E = 1.d0
else
if ( unit .eq. 1 ) then
norm_Th = 1.d0/(EL_CHG*MU_ZERO*1.d20)
else
norm_Th = 1.d0/(K_BOLTZ*MU_ZERO*1.d20)
end if
norm_T = sqrt_mu0_rho0
norm_N = 1.d20
norm_M = 1.d20*MASS_PROTON
norm_E = 1.d0/MU_ZERO
end if
write(*,*) "time:",cor%time
! --- 0 stime
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%time
! --- 1 temperature
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%temperature
do c=1,OUTPUT_MAX
select case( output_commands(c) )
case ('t')
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%time * norm_T
case ('Te')
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%temperature * norm_Th
case ('ne')
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%density * norm_N
case ('fractions')
do iz = 1, cor%n_Z0
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%fractions(iz)
end do
case ('ion_rate')
do iz = 1, cor%ad%n_Z
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%ion_rate(iz) / norm_T / norm_N
end do
case ('rec_rate')
do iz = 1, cor%ad%n_Z
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%rec_rate(iz) / norm_T / norm_N
end do
case ('Z_eff')
z_EFF= cor%z_EFF()
write(pointer,'(a,ES21.15)',advance='no') ' ',z_EFF
case ('P_rad')
P_rad= cor%Prad()
write(pointer,'(a,ES21.15)',advance='no') ' ',P_rad * norm_E / norm_T
case ('L_rad')
L_rad= cor%Lrad()
write(pointer,'(a,ES21.15)',advance='no') ' ',L_rad * norm_E / norm_T / norm_N**2
case ('<Z>')
z_avg= cor%z_avg()
write(pointer,'(a,ES21.15)',advance='no') ' ',z_avg
case ('log10_Te')
log10_Te = log10(cor%temperature/(K_BOLTZ*MU_ZERO*1.d20))
write(pointer,'(a,ES21.15)',advance='no') ' ',log10_Te
case ('log10_ne')
log10_ne = log10(cor%density*1.d20)
write(pointer,'(a,ES21.15)',advance='no') ' ',log10_ne
end select
end do
! --- 2 temperature
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%density
! --- 2+n_Z0 to 3+2*n_Z0-1 f states
do iz = 1, cor%n_Z0
write(pointer,'(a,ES21.15)',advance='no') ' ',cor%fractions(iz)
end do
! --- new line
write(pointer,'(a)') ' '
......@@ -307,6 +399,7 @@ do iz=1,cor%ad%n_Z
call cor%ad%SCD%interp(iz, log10_ne, log10_Te, cor%ion_rate(iz))
call cor%ad%PRB%interp(iz, log10_ne, log10_Te, rad_RB)
call cor%ad%PLT%interp(iz, log10_ne, log10_Te, rad_LT)
cor%rad_RB(iz) = rad_RB ; cor%rad_LT(iz) = rad_LT
cor%rad(iz) = rad_RB + rad_LT
end do
......@@ -314,13 +407,15 @@ end do
cor%rec_rate = cor%rec_rate * sqrt_mu0_rho0 * 1.d20
cor%ion_rate = cor%ion_rate * sqrt_mu0_rho0 * 1.d20
cor%rad = cor%rad * MU_ZERO * sqrt_mu0_rho0 * 1.d40
cor%rad_LT = cor%rad_LT * MU_ZERO * sqrt_mu0_rho0 * 1.d40
cor%rad_RB = cor%rad_RB * MU_ZERO * sqrt_mu0_rho0 * 1.d40
end subroutine calculate_GCRs
! -------------------------------------------------------------------------------------
! --- Calculates Heat source
! --- Calculates Cooling source
! -------------------------------------------------------------------------------------
function calculate_cooling(cor)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
......@@ -362,12 +457,35 @@ end function calculate_equilibrium
! -------------------------------------------------------------------------------------
function calculate_Prad(cor)
class(non_coronal), intent(in) :: cor
real*8 :: calculate_Prad !< Output power in
real*8 :: calculate_Prad !< Output power in W m-3
!calculate_Prad = cor%density*cor%density_imp*dot_product(cor%fractions(2:cor%n_Z0), cor%rad)
calculate_Prad = cor%density*cor%density_imp*(dot_product(cor%fractions(1:cor%ad%n_Z), cor%rad_LT) + dot_product(cor%fractions(2:cor%n_Z0), cor%rad_RB))
calculate_Prad = cor%density*cor%density_imp*dot_product(cor%fractions(2:cor%n_Z0), cor%rad)
end function calculate_Prad
! -------------------------------------------------------------------------------------
! --- Calculates L_rad
! -------------------------------------------------------------------------------------
function calculate_Lrad(cor)
class(non_coronal), intent(in) :: cor
real*8 :: calculate_Lrad !< Output power in W m3
calculate_Lrad = dot_product(cor%fractions(1:cor%ad%n_Z), cor%rad_LT ) + dot_product(cor%fractions(2:cor%n_Z0), cor%rad_RB )
end function calculate_Lrad
! -------------------------------------------------------------------------------------
! --- Calculates average charge state for given f
! -------------------------------------------------------------------------------------
function calculate_z_avg(cor)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
real*8 :: calculate_z_avg
calculate_z_avg = dot_product(cor%fractions,cor%Z)
end function calculate_z_avg
! -------------------------------------------------------------------------------------
......@@ -375,9 +493,12 @@ end function calculate_Prad
! -------------------------------------------------------------------------------------
function calculate_z_EFF(cor)
class(non_coronal), intent(in) :: cor !< Coronal equilibrium type
real*8 :: calculate_z_EFF
real*8 :: calculate_z_EFF, down
calculate_z_EFF = dot_product(cor%fractions,cor%ZZ)
down = dot_product(cor%fractions,cor%Z)
calculate_z_EFF = dot_product(cor%fractions,cor%Z)
if (down .gt. 0.d0 ) calculate_z_EFF = calculate_z_EFF/down
end function calculate_z_EFF
......
......@@ -16,4 +16,6 @@ module constants
real*8, parameter :: MOLE_NUMBER = 6.02214076d23 !< The Avogadro constant
real*8, parameter :: sqrt_mu0_rho0 = 6.483638667465844d-007
integer, parameter :: OUTPUT_MAX = 10
end module constants
File mode changed from 100644 to 100755
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment