Testing charge exchange

Here are the first three steps of charge exchange for a simulation where the densities of stellar and planetary material differ by five orders of magnitude. 

 q(iHIIhot) =    55.0433622655540     
 q(iHcold) =    5504281.18319314     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 After step:
 q(iHIIhot) =    54.7826014982936     
 q(iHcold) =    5478205.36722786     
 q(iHIIcold) =   0.000000000000000E+000
 q(iHhot) =   0.000000000000000E+000
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHhot protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHIIcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Before iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 After iHcold protection:
 q(iHIIhot) =    54.7754817989544     
 q(iHcold) =    5478115.68039737     
 q(iHIIcold) =   6.222884294905237E-003
 q(iHhot) =   6.222884294905237E-003
 Advanced level  0 to tnext= 0.6000E-08 with dt= 0.6000E-08 CFL= 0.5764E-07 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  75%  16%  91% 
 Cell updates/second =        211       209 100%
 Wall Time Remaining =     9.1 day at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2004E+03
 Before step:
 q(iHIIhot) =    54.7752475699149     
 q(iHcold) =    5478092.25506669     
 q(iHIIcold) =   6.222857684820517E-003
 q(iHhot) =   6.222857684820517E-003
 After step:
 q(iHIIhot) =    54.7751651392209     
 q(iHcold) =    5478084.01114326     
 q(iHIIcold) =   6.222848320107427E-003
 q(iHhot) =   6.222848320107427E-003
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHhot protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHIIcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Before iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 After iHcold protection:
 q(iHIIhot) =    22.3818513902974     
 q(iHcold) =    5478051.58934171     
 q(iHIIcold) =    32.3995363123598     
 q(iHhot) =    32.3995363123598     
 Advanced level  0 to tnext= 0.3124E-04 with dt= 0.3124E-04 CFL= 0.3001E-03 max speed= 0.1921E+02 radiation sub-cycles=       1
 Info allocations    =   148.7 kb  148.7 kb
 message allocations =   ------    ------  
 sweep allocations   =   ------    122.5 kb
 filling fractions   = 
 Current efficiency  =  77%  14%  90% 
 Cell updates/second =        303       302 100%
 Wall Time Remaining =     1.8 min at frame    0.0 of     10
 AMR Speed-Up Factor =       0.2879E+03
 Before step:
 q(iHIIhot) =    22.3818513598937     
 q(iHcold) =    5478051.58190027     
 q(iHIIcold) =    32.3995362683480     
 q(iHhot) =    32.3995362683480     

----

 After step:
 q(iHIIhot) =    22.3818513491939     
 q(iHcold) =    5478051.57928144     
 q(iHIIcold) =    32.3995362528591     
 q(iHhot) =    32.3995362528591     
 Charge exchange coefficient is   3.455999999999997E-003
 Before iHIIhot protection:
 q(iHIIhot) =   -4.07465558395285     
 q(iHcold) =    5478025.12276545     
 q(iHIIcold) =    58.8560431859153     
 q(iHhot) =    58.8560431859153</span>

----

 After iHIIhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     

----

 Before iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHhot protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHIIcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 Before iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624     
 After iHcold protection:
 q(iHIIhot) =   0.000000000000000E+000
 q(iHcold) =    5478029.19742104     
 q(iHIIcold) =    54.7813876019624     
 q(iHhot) =    54.7813876019624
 Advanced level  0 to tnext= 0.9368E-04 with dt= 0.6244E-04...

Jonathan was right on the money with the timestepping being too long. I've broken out the relevant blocks above. For the charge exchange step that follows the first line, we get an overall of

Since HIIhot is only ~22, this gives us HIIhot of ~-4 after charge exchange is complete.

This run is from after I implemented protections for any species less than 0. The section that follows the first line shows that 4 is added back to HIIhot and Hcold, and removed from HIIcold and Hhot, reversing the excess charge exchange that brought HIIhot negative.

This method of protection can result in large oscillations arount equilibrium. Better would be to prevent any one species from changing more than e.g. 10%, as we do with the other quantities, but this is more computationally expensive thanks to the amount of subcycling that would be required in some places (namely those with large in comparison to any one species).

Comments

No comments.