From 02a0624b32fa54c3091633374e4ae70e8e2a1b8a Mon Sep 17 00:00:00 2001 From: Alexander Nozik Date: Thu, 6 Dec 2018 17:32:01 +0300 Subject: [PATCH] Separate generator for each track --- .../inr/numass/trapping/SimulationManager.kt | 29 ++++++---- .../kotlin/inr/numass/trapping/Simulator.kt | 56 ++++++++----------- .../kotlin/inr/numass/trapping/Trapping.kt | 5 +- 3 files changed, 43 insertions(+), 47 deletions(-) diff --git a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt index b90d5b9..8a62bf1 100644 --- a/src/main/kotlin/inr/numass/trapping/SimulationManager.kt +++ b/src/main/kotlin/inr/numass/trapping/SimulationManager.kt @@ -1,14 +1,17 @@ package inr.numass.trapping 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.simple.RandomSource import java.io.File import java.io.PrintStream import java.nio.file.Files 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. @@ -20,7 +23,10 @@ class SimulationManager() { 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 initialE = 18000.0 @@ -85,8 +91,8 @@ class SimulationManager() { 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() @@ -100,18 +106,19 @@ class SimulationManager() { val outputPath = outputDirectory.toPath().resolve("$fileName.out") + PrintStream(Files.newOutputStream(outputPath, StandardOpenOption.CREATE, StandardOpenOption.TRUNCATE_EXISTING, StandardOpenOption.WRITE)).use { output -> 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()) - output.printf("%s\t%s\t%s\t%s\t%s\t%s%n", "E", "theta", "theta_start", "colNum", "L", "state") - Stream - .generate { + output.printf("%s\t%s\t%s\t%s\t%s\t%s\t%s%n", "id", "E", "theta", "theta_start", "colNum", "L", "state") + LongStream.rangeClosed(1, num.toLong()).parallel() + .mapToObj { id -> + val generator = RandomGeneratorBridge(generatorFactory(id)) val theta = Math.acos(1 - 2 * generator.nextDouble())// from 0 to Pi 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) .forEach { res -> printOne(output, res) @@ -148,7 +155,7 @@ class SimulationManager() { } 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() } diff --git a/src/main/kotlin/inr/numass/trapping/Simulator.kt b/src/main/kotlin/inr/numass/trapping/Simulator.kt index 90f135c..6ce096b 100644 --- a/src/main/kotlin/inr/numass/trapping/Simulator.kt +++ b/src/main/kotlin/inr/numass/trapping/Simulator.kt @@ -23,8 +23,8 @@ class Simulator( val thetaPinch: Double, val gasDensity: Double,// m^-3 val bSource: Double, - val magneticField: ((Double) -> Double)?, - val generator: RandomGenerator) { + val magneticField: ((Double) -> Double)? +) { enum class EndState { @@ -55,20 +55,12 @@ class Simulator( //проверяем нормировку assert(Precision.equals(sigmaEl + sigmaExc + sigmaIon, 1.0, 1e-2)) - val alpha = generator.nextDouble() + val alpha = pos.generator.nextDouble() - val delta: Pair - if (alpha > sigmaEl) { - if (alpha > sigmaEl + sigmaExc) { - //ionization case - delta = Scatter.randomion(pos.e, generator) - } else { - //excitation case - delta = Scatter.randomexc(pos.e, generator) - } - } else { - // elastic - delta = Scatter.randomel(pos.e, generator) + val delta: Pair = when { + alpha < sigmaEl -> Scatter.randomel(pos.e, pos.generator)// elastic + alpha > sigmaEl + sigmaExc -> Scatter.randomion(pos.e, pos.generator)//ionization case + else -> Scatter.randomexc(pos.e, pos.generator)//excitation case } //Обновляем значени угла и энергии независимо ни от чего @@ -85,7 +77,7 @@ class Simulator( * @param e * @return */ - private fun freePath(e: Double): Double { + private fun freePath(generator: RandomGenerator, e: Double): Double { //FIXME redundant cross-section calculation //All cross-sections are in m^2 return ExponentialDistribution(generator, 1.0 / Scatter.sigmaTotal(e) / gasDensity).sample() @@ -196,15 +188,15 @@ class Simulator( * Симулируем один пробег электрона от начального значения и до вылетания из * иточника или до того момента, как энергия становится меньше 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(initTheta > 0 && initTheta < Math.PI) assert(abs(initZ) <= SOURCE_LENGTH / 2.0) - val pos = State(initEnergy, initTheta, initZ) + val pos = State(generator, initEnergy, initTheta, initZ) 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(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() { @@ -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! @@ -245,7 +237,7 @@ class Simulator( * @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 */ - 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 */ @@ -349,7 +341,7 @@ class Simulator( /** * Reverse electron direction */ - internal fun flip() { + fun flip() { if (theta < 0 || theta > PI) { throw Error() } @@ -361,7 +353,7 @@ class Simulator( * * @return */ - internal fun field(): Double { + fun field(): Double { return this@Simulator.field(z) } @@ -370,7 +362,7 @@ class Simulator( * * @return */ - internal fun realTheta(): Double { + fun realTheta(): Double { if (magneticField == null) { return theta } else { @@ -391,7 +383,7 @@ class Simulator( * @param dTheta * @return resulting angle */ - internal fun addTheta(dTheta: Double): Double { + fun addTheta(dTheta: Double): Double { //Генерируем случайный фи val phi = generator.nextDouble() * 2.0 * Math.PI @@ -406,7 +398,7 @@ class Simulator( val result = rot.applyTo(init.cartesian) - val newtheta = acos(result.z) + val newTheta = acos(result.z) // //следим чтобы угол был от 0 до Pi // if (newtheta < 0) { @@ -418,10 +410,10 @@ class Simulator( //change back to virtual angles if (magneticField == null) { - theta = newtheta + theta = newTheta } else { - theta = asin(sin(newtheta) * sqrt(bSource / field())) - if (newtheta > PI / 2) { + theta = asin(sin(newTheta) * sqrt(bSource / field())) + if (newTheta > PI / 2) { theta = PI - theta } } @@ -435,8 +427,8 @@ class Simulator( companion object { - val SOURCE_LENGTH = 3.0 - private val DELTA_L = 0.1 //step for propagate calculation + const val SOURCE_LENGTH = 3.0 + private const val DELTA_L = 0.1 //step for propagate calculation } diff --git a/src/main/kotlin/inr/numass/trapping/Trapping.kt b/src/main/kotlin/inr/numass/trapping/Trapping.kt index eb1a61f..9ed354d 100644 --- a/src/main/kotlin/inr/numass/trapping/Trapping.kt +++ b/src/main/kotlin/inr/numass/trapping/Trapping.kt @@ -1,7 +1,5 @@ package inr.numass.trapping -import org.apache.commons.rng.simple.RandomSource -import java.io.File import java.time.Duration import java.time.Instant @@ -21,12 +19,11 @@ fun main(args: Array) { SimulationManager().apply { comment = "Out of the box cross-sections" fileName = "trap[$e]" - generator = RandomSource.create(RandomSource.MWC_256) setFields(0.6, 3.7, 7.2) gasDensity = 1e19 initialE = e range = 4000.0 - }.simulateAll(1e7) + }.simulateAll(1_000_000) } val finishTime = Instant.now()