NPM/trapping
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Separate generator for each track

Browse Source
This commit is contained in:
Alexander Nozik 2018-12-06 17:32:01 +03:00
parent 839526c7af
commit 02a0624b32
3 changed files with 43 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
src/main/kotlin/inr/numass/trapping/SimulationManager.kt
View File

@ -1,14 +1,17 @@
package inr.numass.trapping package inr.numass.trapping
import org.apache.commons.math3.analysis.interpolation.LinearInterpolator import org.apache.commons.math3.analysis.interpolation.LinearInterpolator
import org.apache.commons.math3.random.SynchronizedRandomGenerator import org.apache.commons.math3.random.JDKRandomGenerator
import org.apache.commons.rng.UniformRandomProvider import org.apache.commons.rng.UniformRandomProvider
import org.apache.commons.rng.simple.RandomSource import org.apache.commons.rng.simple.RandomSource
import java.io.File import java.io.File
import java.io.PrintStream import java.io.PrintStream
import java.nio.file.Files import java.nio.file.Files
import java.nio.file.StandardOpenOption import java.nio.file.StandardOpenOption
import java.util.stream.Stream import java.util.stream.LongStream
private val seedGenerator = JDKRandomGenerator()
/** /**
* Created by darksnake on 04-Jun-16. * Created by darksnake on 04-Jun-16.
@ -20,7 +23,10 @@ class SimulationManager() {
var comment = "" var comment = ""
var generator: UniformRandomProvider = RandomSource.create(RandomSource.SPLIT_MIX_64) /**
* A supplier for random generator. Each track has its own generator
*/
var generatorFactory: (Long) -> UniformRandomProvider = { RandomSource.create(RandomSource.MT_64, seedGenerator.nextInt()) }
var reportFilter: (Simulator.SimulationResult) -> Boolean = { it.state == Simulator.EndState.ACCEPTED } var reportFilter: (Simulator.SimulationResult) -> Boolean = { it.state == Simulator.EndState.ACCEPTED }
var initialE = 18000.0 var initialE = 18000.0
@ -85,8 +91,8 @@ class SimulationManager() {
outputDirectory.mkdirs() outputDirectory.mkdirs()
} }
val generator = SynchronizedRandomGenerator(RandomGeneratorBridge(generator))
val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField, generator) val simulator = Simulator(eLow, thetaTransport, thetaPinch, gasDensity, bSource, magneticField)
val counter = Counter() val counter = Counter()
@ -100,18 +106,19 @@ class SimulationManager() {
val outputPath = outputDirectory.toPath().resolve("$fileName.out") val outputPath = outputDirectory.toPath().resolve("$fileName.out")
PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output -> PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output ->
output.println("# " + header.replace("\n", "\n# "))//adding comment symbols output.println("# " + header.replace("\n", "\n# "))//adding comment symbols
System.out.printf("%nStarting simulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong()) System.out.printf("%nStarting simulation with initial energy %g and %d electrons.%n%n", initialE, num.toLong())
output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state") output.printf("%s\t%s\t%s\t%s\t%s\t%s\t%s%n", "id", "E", "theta", "theta_start", "colNum", "L", "state")
Stream LongStream.rangeClosed(1, num.toLong()).parallel()
.generate { .mapToObj { id ->
val generator = RandomGeneratorBridge(generatorFactory(id))
val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi
val z = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH val z = (generator.nextDouble() - 0.5) * Simulator.SOURCE_LENGTH
simulator.simulate(initialE, theta, z).also { counter.count(it) } simulator.simulate(id, generator, initialE, theta, z).also { counter.count(it) }
} }
.limit(num.toLong()).parallel()
.filter(reportFilter) .filter(reportFilter)
.forEach { res -> .forEach { res ->
printOne(output, res) printOne(output, res)
@ -148,7 +155,7 @@ class SimulationManager() {
} }
private fun printOne(out: PrintStream, res: Simulator.SimulationResult) { private fun printOne(out: PrintStream, res: Simulator.SimulationResult) {
out.printf("%g\t%g\t%g\t%d\t%g\t%s%n", res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString()) out.printf("%d\t%g\t%g\t%g\t%d\t%g\t%s%n", res.id, res.E, res.theta * 180 / Math.PI, res.initTheta * 180 / Math.PI, res.collisionNumber, res.l, res.state.toString())
out.flush() out.flush()
} }

56
src/main/kotlin/inr/numass/trapping/Simulator.kt
View File

@ -23,8 +23,8 @@ class Simulator(
val thetaPinch: Double, val thetaPinch: Double,
val gasDensity: Double,// m^-3 val gasDensity: Double,// m^-3
val bSource: Double, val bSource: Double,
val magneticField: ((Double) -> Double)?, val magneticField: ((Double) -> Double)?
val generator: RandomGenerator) { ) {
enum class EndState { enum class EndState {
@ -55,20 +55,12 @@ class Simulator(
//проверяем нормировку //проверяем нормировку
assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2)) assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2))
val alpha = generator.nextDouble() val alpha = pos.generator.nextDouble()
val delta: Pair<Double, Double> val delta: Pair<Double, Double> = when {
if (alpha > sigmaEl) { alpha < sigmaEl -> Scatter.randomel(pos.e, pos.generator)// elastic
if (alpha > sigmaEl + sigmaExc) { alpha > sigmaEl + sigmaExc -> Scatter.randomion(pos.e, pos.generator)//ionization case
//ionization case else -> Scatter.randomexc(pos.e, pos.generator)//excitation case
delta = Scatter.randomion(pos.e, generator)
} else {
//excitation case
delta = Scatter.randomexc(pos.e, generator)
}
} else {
// elastic
delta = Scatter.randomel(pos.e, generator)
} }
//Обновляем значени угла и энергии независимо ни от чего //Обновляем значени угла и энергии независимо ни от чего
@ -85,7 +77,7 @@ class Simulator(
* @param e * @param e
* @return * @return
*/ */
private fun freePath(e: Double): Double { private fun freePath(generator: RandomGenerator, e: Double): Double {
//FIXME redundant cross-section calculation //FIXME redundant cross-section calculation
//All cross-sections are in m^2 //All cross-sections are in m^2
return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample() return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample()
@ -196,15 +188,15 @@ class Simulator(
* Симулируем один пробег электрона от начального значения и до вылетания из * Симулируем один пробег электрона от начального значения и до вылетания из
* иточника или до того момента, как энергия становится меньше eLow. * иточника или до того момента, как энергия становится меньше eLow.
*/ */
fun simulate(initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult { fun simulate(id: Long, generator: RandomGenerator, initEnergy: Double, initTheta: Double, initZ: Double): SimulationResult {
assert(initEnergy > 0) assert(initEnergy > 0)
assert(initTheta > 0 && initTheta < Math.PI) assert(initTheta > 0 && initTheta < Math.PI)
assert(abs(initZ) <= SOURCE_LENGTH / 2.0) assert(abs(initZ) <= SOURCE_LENGTH / 2.0)
val pos = State(initEnergy, initTheta, initZ) val pos = State(generator, initEnergy, initTheta, initZ)
while (!pos.isFinished) { while (!pos.isFinished) {
val dl = freePath(pos.e) // path to next scattering val dl = freePath(generator, pos.e) // path to next scattering
// propagate to next scattering position // propagate to next scattering position
propagate(pos, dl) propagate(pos, dl)
@ -220,7 +212,7 @@ class Simulator(
} }
} }
return SimulationResult(pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l) return SimulationResult(id, pos.endState, pos.e, pos.theta, initTheta, pos.colNum, pos.l)
} }
fun resetDebugCounters() { fun resetDebugCounters() {
@ -236,7 +228,7 @@ class Simulator(
} }
} }
class SimulationResult(var state: EndState, var E: Double, var theta: Double, var initTheta: Double, var collisionNumber: Int, var l: Double) data class SimulationResult(val id: Long, val state: EndState, val E: Double, val theta: Double, val initTheta: Double, val collisionNumber: Int, val l: Double)
/** /**
* Current electron position in simulation. Not thread safe! * Current electron position in simulation. Not thread safe!
@ -245,7 +237,7 @@ class Simulator(
* @property theta Current theta recalculated to the field in the center of the source * @property theta Current theta recalculated to the field in the center of the source
* @property z current z. Zero is the center of the source * @property z current z. Zero is the center of the source
*/ */
private inner class State(internal var e: Double, internal var theta: Double, internal var z: Double) { private inner class State(val generator: RandomGenerator, internal var e: Double, internal var theta: Double, internal var z: Double) {
/** /**
* Current total path * Current total path
*/ */
@ -349,7 +341,7 @@ class Simulator(
/** /**
* Reverse electron direction * Reverse electron direction
*/ */
internal fun flip() { fun flip() {
if (theta < 0 || theta > PI) { if (theta < 0 || theta > PI) {
throw Error() throw Error()
} }
@ -361,7 +353,7 @@ class Simulator(
* *
* @return * @return
*/ */
internal fun field(): Double { fun field(): Double {
return this@Simulator.field(z) return this@Simulator.field(z)
} }
@ -370,7 +362,7 @@ class Simulator(
* *
* @return * @return
*/ */
internal fun realTheta(): Double { fun realTheta(): Double {
if (magneticField == null) { if (magneticField == null) {
return theta return theta
} else { } else {
@ -391,7 +383,7 @@ class Simulator(
* @param dTheta * @param dTheta
* @return resulting angle * @return resulting angle
*/ */
internal fun addTheta(dTheta: Double): Double { fun addTheta(dTheta: Double): Double {
//Генерируем случайный фи //Генерируем случайный фи
val phi = generator.nextDouble() * 2.0 * Math.PI val phi = generator.nextDouble() * 2.0 * Math.PI
@ -406,7 +398,7 @@ class Simulator(
val result = rot.applyTo(init.cartesian) val result = rot.applyTo(init.cartesian)
val newtheta = acos(result.z) val newTheta = acos(result.z)
// //следим чтобы угол был от 0 до Pi // //следим чтобы угол был от 0 до Pi
// if (newtheta < 0) { // if (newtheta < 0) {
@ -418,10 +410,10 @@ class Simulator(
//change back to virtual angles //change back to virtual angles
if (magneticField == null) { if (magneticField == null) {
theta = newtheta theta = newTheta
} else { } else {
theta = asin(sin(newtheta) * sqrt(bSource / field())) theta = asin(sin(newTheta) * sqrt(bSource / field()))
if (newtheta > PI / 2) { if (newTheta > PI / 2) {
theta = PI - theta theta = PI - theta
} }
} }
@ -435,8 +427,8 @@ class Simulator(
companion object { companion object {
val SOURCE_LENGTH = 3.0 const val SOURCE_LENGTH = 3.0
private val DELTA_L = 0.1 //step for propagate calculation private const val DELTA_L = 0.1 //step for propagate calculation
} }

5
src/main/kotlin/inr/numass/trapping/Trapping.kt
View File

@ -1,7 +1,5 @@
package inr.numass.trapping package inr.numass.trapping
import org.apache.commons.rng.simple.RandomSource
import java.io.File
import java.time.Duration import java.time.Duration
import java.time.Instant import java.time.Instant
@ -21,12 +19,11 @@ fun main(args: Array<String>) {
SimulationManager().apply { SimulationManager().apply {
comment = "Out of the box cross-sections" comment = "Out of the box cross-sections"
fileName = "trap[$e]" fileName = "trap[$e]"
generator = RandomSource.create(RandomSource.MWC_256)
setFields(0.6, 3.7, 7.2) setFields(0.6, 3.7, 7.2)
gasDensity = 1e19 gasDensity = 1e19
initialE = e initialE = e
range = 4000.0 range = 4000.0
}.simulateAll(1e7) }.simulateAll(1_000_000)
} }
val finishTime = Instant.now() val finishTime = Instant.now()