diff options
| -rw-r--r-- | xmds2/realistic_Rb/Makefile | 3 | ||||
| -rwxr-xr-x | xmds2/realistic_Rb/tests/run_tests.py | 350 | ||||
| -rw-r--r-- | xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds | 667 | ||||
| -rw-r--r-- | xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil | 699 | ||||
| -rw-r--r-- | xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat | bin | 0 -> 3240832 bytes |
5 files changed, 1719 insertions, 0 deletions
diff --git a/xmds2/realistic_Rb/Makefile b/xmds2/realistic_Rb/Makefile index e088e15..34d73e1 100644 --- a/xmds2/realistic_Rb/Makefile +++ b/xmds2/realistic_Rb/Makefile @@ -46,6 +46,9 @@ $(png_targets): %.png : %.pdf pdf: $(pdf_targets) +test: + cd tests; ./run_tests.py + $(pdf_targets): %.pdf : %.eps cat $< | ps2eps -B > __tt.eps epspdf __tt.eps $@ diff --git a/xmds2/realistic_Rb/tests/run_tests.py b/xmds2/realistic_Rb/tests/run_tests.py new file mode 100755 index 0000000..61f4eea --- /dev/null +++ b/xmds2/realistic_Rb/tests/run_tests.py @@ -0,0 +1,350 @@ +#!/usr/bin/env python +# encoding: utf-8 +""" +run_tests.py + +Created by Graham Dennis on 2008-06-15. + +Copyright (c) 2008-2012, Graham Dennis + +This program is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program. If not, see <http://www.gnu.org/licenses/>. + +""" + +import xpdeint.Python24Support + +import os +import re +import sys +import getopt +import shutil +import hashlib +import unittest +import subprocess + +from xml.dom import minidom +import xpdeint.minidom_extras +from xpdeint import CodeParser + +from xpdeint.XSILFile import XSILFile + +import numpy + +help_message = ''' +The help message goes here. +''' + + +class Usage(Exception): + def __init__(self, msg): + self.msg = msg + +def pass_nan_test(array1, array2): + """Return `True` if isNaN(`array1`) == isNaN(`array2`)""" + # NaN test. array2 is allowed to be NaN at an index if array1 is also NaN there. + nanTestPassed = numpy.equal(numpy.isnan(array1), numpy.isnan(array2)).all() + return nanTestPassed + +def array_approx_equal(array1, array2, absTol, relTol): + """Return `True` if all of (`array1` - `array2`) <= `absTol` or (`array1` - `array2`) <= `relTol` * `array2`""" + diff = array1-array2 + # NaN values would fail this test. So we have to exclude them. But only exclude them if array2 (the expected results) + # have a NaN + return numpy.logical_or(numpy.logical_or(numpy.abs(diff) <= 0.5 * relTol * (numpy.abs(array2) + numpy.abs(array1)), numpy.abs(diff) <= absTol), numpy.isnan(array2)).all() + +def scriptTestingFunction(root, scriptName, testDir, absPath, self): + if not os.path.exists(testDir): + os.makedirs(testDir) + + proc = subprocess.Popen('xmds2 --no-version ' + absPath, + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + cwd=testDir) + (stdout, stderr) = proc.communicate() + returnCode = proc.wait() + + message = ''.join(["\n%(handleName)s:\n%(content)s" % locals() for handleName, content in [('stdout', stdout), ('stderr', stderr)] if content]) + + # A few tests require XMDS1. If XMDS1 isn't present we should just + # skip that test rather than failing. + # The skip functionality for the unittest class is only available + # in python 2.7 and later, so check for that too. + if returnCode != 0 and sys.version_info[:2] >= (2, 7): + if re.search(r'^The missing \w+ feature\(s\) were: .*xmds.*', message, re.MULTILINE): + self.skipTest("Skipping test as XMDS1 is required and not installed") + + # A few tests require specific features. If it isn't available, skip the test + # rather than failing. + # The skip functionality for the unittest class is only available + # in python 2.7 and later, so check for that too. + if returnCode != 0 and sys.version_info[:2] >= (2, 7): + if re.search(r'^The missing \w+ feature\(s\) were:', message, re.MULTILINE): + self.skipTest("Skipping test as feature required is not installed") + + self.assert_(returnCode == 0, ("Failed to compile." % locals()) + message) + + xmlDocument = minidom.parse(absPath) + simulationElement = xmlDocument.getChildElementByTagName('simulation') + nameElement = simulationElement.getChildElementByTagName('name') + testingElement = simulationElement.getChildElementByTagName('testing') + + simulationName = nameElement.innerText() + + # If the source is the same as the last known good, then we don't need to compile or execute the simulation. + sourceFilePath = os.path.join(testDir, simulationName + '.cc') + checksumFilePath = os.path.join(testDir, simulationName + '_last_known_good.checksum') + sourceContents = file(sourceFilePath).read() + h = hashlib.sha1() + h.update(sourceContents) + currentChecksum = h.hexdigest() + + if os.path.exists(checksumFilePath): + lastKnownGoodChecksum = file(checksumFilePath).read() + + if lastKnownGoodChecksum == currentChecksum: + # The checksums check out, so we don't need to go any further + return + + # Now we have compiled, we need to copy any input data needed and then run the simulation + inputXSILElements = testingElement.getChildElementsByTagName('input_xsil_file', optional=True) + + filesToCopy = [] + + for inputXSILElement in inputXSILElements: + name = inputXSILElement.getAttribute('name').strip() + filesToCopy.append(name) + inputXSILFile = XSILFile(os.path.join(os.path.split(absPath)[0], name), loadData=False) + filesToCopy.extend([os.path.join(os.path.split(name)[0], xsil.data.filename) for xsil in inputXSILFile.xsilObjects if hasattr(xsil.data, 'filename')]) + + for fileToCopy in filesToCopy: + sourceFile = os.path.join(os.path.split(absPath)[0], fileToCopy) + shutil.copy(sourceFile, testDir) + + # Allow command-line arguments to be specified for the simulation + commandLineElement = testingElement.getChildElementByTagName('command_line', optional=True) + argumentsElement = testingElement.getChildElementByTagName('arguments', optional=True) + commandLineString = './' + simulationName + if commandLineElement: + # The command line element overrides the prefix + commandLineString = commandLineElement.innerText().strip() + if argumentsElement: + commandLineString += ' ' + argumentsElement.innerText().strip() + + simulationProc = subprocess.Popen(commandLineString, + shell=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + cwd=testDir) + (stdout, stderr) = simulationProc.communicate() + returnCode = simulationProc.wait() + + self.assert_(returnCode == 0, "Failed to execute compiled simulation correctly." % locals()) + + # The next thing to check is that the generated data agrees with the expected data to within the set error margins. + xsilFileElements = testingElement.getChildElementsByTagName('xsil_file', optional=True) + for xsilFileElement in xsilFileElements: + sourceFile = xsilFileElement.getAttribute('name').strip() + expectedResultsFile = xsilFileElement.getAttribute('expected').strip() + # Defaults + absoluteTolerance = 0 + relativeTolerance = 1e-9 + + if xsilFileElement.hasAttribute('absolute_tolerance'): + absoluteTolerance = float(xsilFileElement.getAttribute('absolute_tolerance')) + if xsilFileElement.hasAttribute('relative_tolerance'): + relativeTolerance = float(xsilFileElement.getAttribute('relative_tolerance')) + + resultsFullPath = os.path.join(testDir, sourceFile) + results = XSILFile(resultsFullPath) + expectedResultsFullPath = os.path.join(os.path.split(absPath)[0], expectedResultsFile) + if not os.path.exists(expectedResultsFullPath): + print >> sys.stderr, "Expected results file '%(expectedResultsFile)s' missing. Using current. " % locals() + + # If there are any NaN's in the results, issue a warning. + for mgNum, o in enumerate(results.xsilObjects): + for v in o.independentVariables: + if numpy.isnan(v['array']).any(): + print >> sys.stderr, "Warning: Coordinate '%s' in moment group %i of file '%s' contains a NaN." % (v['name'], mgNum+1, sourceFile) + for v in o.dependentVariables: + if numpy.isnan(v['array']).any(): + print >> sys.stderr, "Warning: Dependent variable '%s' in moment group %i of file '%s' contains a NaN." % (v['name'], mgNum+1, sourceFile) + + resultsFileContents = file(resultsFullPath).read() + + for xsilObject in results.xsilObjects: + if hasattr(xsilObject.data, 'filename'): + # If the moment group has a data file name, then we need to copy it to the expected results file + newDataFilename = xsilObject.data.filename.replace(os.path.splitext(sourceFile)[0], os.path.splitext(expectedResultsFile)[0], 1) + + resultsFileContents = resultsFileContents.replace(xsilObject.data.filename, newDataFilename) + + shutil.copyfile(os.path.join(testDir, xsilObject.data.filename), + os.path.join(os.path.split(absPath)[0], newDataFilename)) + + file(expectedResultsFullPath, 'w').write(resultsFileContents) + else: + expectedResults = XSILFile(expectedResultsFullPath) + + self.assert_(len(results.xsilObjects) == len(expectedResults.xsilObjects)) + + momentGroupElements = xsilFileElement.getChildElementsByTagName('moment_group', optional=True) + if momentGroupElements: + self.assert_(len(momentGroupElements) == len(results.xsilObjects)) + else: + momentGroupElements = [None]*len(results.xsilObjects) + + for mgNum, (o1, o2, mgElem) in enumerate(zip(results.xsilObjects, expectedResults.xsilObjects, momentGroupElements)): + currentAbsoluteTolerance = absoluteTolerance + currentRelativeTolerance = relativeTolerance + self.assert_(len(o1.independentVariables) == len(o2.independentVariables), + "The number of independent variables in moment group %(mgNum)i doesn't match." % locals()) + self.assert_(len(o1.dependentVariables) == len(o2.dependentVariables), + "The number of dependent variables in moment group %(mgNum)i doesn't match." % locals()) + + if mgElem: + if mgElem.hasAttribute('absolute_tolerance'): + currentAbsoluteTolerance = float(mgElem.getAttribute('absolute_tolerance')) + if mgElem.hasAttribute('relative_tolerance'): + currentRelativeTolerance = float(mgElem.getAttribute('relative_tolerance')) + + self.assert_(currentAbsoluteTolerance != None and currentRelativeTolerance != None, "An absolute and a relative tolerance must be specified.") + + for v1, v2 in zip(o1.independentVariables, o2.independentVariables): + self.assert_(v1['name'] == v2['name']) + self.assert_(v1['length'] == v2['length']) + # These are the coordinates, we just specify a constant absolute and relative tolerance. + # No-one should need to change these + self.assert_(array_approx_equal(v1['array'], v2['array'], 1e-7, 1e-6), + "Coordinate '%s' in moment group %i of file '%s' didn't pass tolerance criteria." % (v1['name'], mgNum+1, sourceFile)) + + for v1, v2 in zip(o1.dependentVariables, o2.dependentVariables): + self.assert_(v1['name'] == v2['name']) + self.assert_(pass_nan_test(v1['array'], v2['array']), + "Dependent variable '%s' in moment group %i of file '%s' had a NaN where the expected results didn't (or vice-versa)." % (v1['name'], mgNum+1, sourceFile)) + self.assert_(array_approx_equal(v1['array'], v2['array'], currentAbsoluteTolerance, currentRelativeTolerance), + "Dependent variable '%s' in moment group %i of file '%s' failed to pass tolerance criteria." % (v1['name'], mgNum+1, sourceFile)) + + # Test has succeeded, so save our checksum for the source file and copy the source file + file(checksumFilePath, 'w').write(currentChecksum) + + lastKnownGoodSourcePath = os.path.join(testDir, simulationName + '_last_known_good.cc') + file(lastKnownGoodSourcePath, 'w').write(sourceContents) + +def partial(func, *args, **keywords): + def newfunc(*fargs, **fkeywords): + newkeywords = keywords.copy() + newkeywords.update(fkeywords) + return func(*(args + fargs), **newkeywords) + return newfunc + + +def main(argv=None): + if argv is None: + argv = sys.argv + try: + try: + opts, args = getopt.getopt(argv[1:], "ho:v", ["help", "output="]) + except getopt.error, msg: + raise Usage(msg) + + # option processing + for option, value in opts: + if option == "-v": + verbose = True + if option in ("-h", "--help"): + raise Usage(help_message) + if option in ("-o", "--output"): + output = value + + except Usage, err: + print >> sys.stderr, sys.argv[0].split("/")[-1] + ": " + str(err.msg) + print >> sys.stderr, "\t for help use --help" + return 2 + + basePath = os.path.dirname(__file__) + + resultsPath = os.path.join(basePath, 'testsuite_results') + if not os.path.exists(resultsPath): + os.mkdir(resultsPath) + resultsPath = os.path.abspath(resultsPath) + + print "Saving test results in %(resultsPath)s" % locals() + + testsuites = {} + baseSuiteName = 'testsuite' + baseSuitePath = os.path.join(basePath, baseSuiteName) + + for root, dirs, files in os.walk(baseSuitePath): + # First remove directories we don't want to traverse + for dirName in ['.svn']: + if dirName in dirs: + dirs.remove(dirName) + # Remove the 'testsuite/' part of the path + dirRelativeToBase = root[(len(baseSuitePath)+1):] + if dirRelativeToBase: + testSuiteName = os.path.join(baseSuiteName, dirRelativeToBase) + else: + testSuiteName = baseSuiteName + + # If we have .xmds files in this path, then create a TestCase subclass + xmdsTestScripts = [filename for filename in files if os.path.splitext(filename)[1].lower() == '.xmds'] + + if xmdsTestScripts: + class ScriptTestCase(unittest.TestCase): + # Create test functions for each test script using 'scriptTestingFunction' + # These test function names are of the form 'test_ScriptName' + for scriptName in xmdsTestScripts: + prefix = os.path.splitext(scriptName)[0] + absPath = os.path.abspath(os.path.join(root, scriptName)) + testDir = os.path.join(resultsPath, dirRelativeToBase) + locals()['test_' + prefix] = partial(scriptTestingFunction, root, scriptName, testDir, absPath) + locals()['test_' + prefix].__doc__ = os.path.join(dirRelativeToBase, scriptName) + + # Create a TestSuite from that class + suite = unittest.defaultTestLoader.loadTestsFromTestCase(ScriptTestCase) + testsuites[testSuiteName] = suite + + if not testSuiteName in testsuites: + testsuites[testSuiteName] = unittest.TestSuite() + + suite = testsuites[testSuiteName] + # Add our TestSuite as a sub-suite of all parent suites + head = testSuiteName + while True: + head, tail = os.path.split(head) + if not head or not tail: + break + testsuites[head].addTest(suite) + + + suitesToRun = list() + if len(args): + for suiteName in args: + fullSuiteName = os.path.join(baseSuiteName, suiteName) + if fullSuiteName in testsuites: + suitesToRun.append(testsuites[fullSuiteName]) + else: + print >> sys.stderr, "Unable to find test '%(suiteName)s'" % locals() + else: + suitesToRun.append(testsuites[baseSuiteName]) + suitesToRun.append(unittest.defaultTestLoader.loadTestsFromModule(CodeParser)) + + fullSuite = unittest.TestSuite(tests=suitesToRun) + + return not unittest.TextTestRunner().run(fullSuite).wasSuccessful() + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds new file mode 100644 index 0000000..8963ae5 --- /dev/null +++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb.xmds @@ -0,0 +1,667 @@ +<?xml version="1.0"?> +<simulation xmds-version="2"> + <testing> + <arguments> --Ndens=1e15 --Lcell=10.0e-2 --Temperature=1e-3 --Pwidth=0.4e-6 --delta1=0 --delta2=0 --delta3=0 --E1o=0 --E2o=1e2 --E3o=0 --E4o=0</arguments> + <xsil_file name="realistic_Rb.xsil" expected="realistic_Rb_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5"> + <moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" /> + </xsil_file> + </testing> + + <name>realistic_Rb</name> + + <author>Eugeniy Mikhailov</author> + <description> + License GPL. + + Solving simplified Rb atom model + with fields propagation along spatial axis Z + with Doppler broadening. + + + We assume four-wave mixing condition when w3-w4=w2-w1 i.e. fields E3 and E4 drive the same + resonance as fields E2 and E1. + + + * --------------- | F=1, 2P_3/2 > + * \ \ + * \ E3_r \ -------- | F=2, 2P_+1/2 > + * \ E4_r \ / \ + * \ \ / E2_l \ + * \ / \ E1_l + * | F=2, 2S_1/2 > -------------- \ + * \ \ + * \ \ + * ------------- | F=1, 2S_1/2 > + * + + + We are solving + dE/dz+(1/c)*dE/dt=i*eta*rho_ij, where j level is higher then i. + Note that E is actually a Rabi frequency of electromagnetic field not the EM field + in xmds terms it looks like + dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space + + VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to + normalize it to something else look drho/dt equation. + No need to renormalizes eta as long as its express through + the upper level decay rate in the same units as Rabi frequency. + </description> + + <features> + <globals> + <![CDATA[ + // Some numerical constants + const double pi = M_PI; + // proportional to splitting ratios sqrt(6) , sqrt(3), sqrt(2) + const double rt6 = 2.449489742783178; + const double rt3 = 1.7320508075688772; + const double rt2 = 1.4142135623730951; + + + const double c=3.e8; + const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K] + const double lambda=794.7e-9; //wavelength in m + // Fields k-vector + const double Kvec = 2*M_PI/lambda; + // Simplified k-vectors + const double Kvec1 = Kvec, Kvec2=Kvec, Kvec3=Kvec; + + const double Gamma_super=6*(2*M_PI*1e6); // characteristic decay rate of upper level used for eta calculations expressed in [1/s] + // eta will be calculated in the <arguments> section + double eta = 0; // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)] + double eta1=0, eta2=0, eta3=0; + + // --------- Atom and cell properties ------------------------- + // range of Maxwell distribution atomic velocities + const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] + // above mass expression is written as (expression is isotopic_mass * atomic_mass_unit) + + // Average sqrt(v^2) in Maxwell distribution for one dimension + // Maxwell related parameters will be calculated in <arguments> section + double v_thermal_averaged=0; + // Maxwell distribution velocities range to take in account in [m/s] + double V_maxwell_min = 0, V_maxwell_max = 0; + + // repopulation rate (atoms flying in/out the laser beam) in [1/s] + const double gt=0.01 *(2*M_PI*1e6); + + // Natural linewidth of j's level in [1/s] + const double g1 = 3.612847284945266e7; + const double g2 = 3.8117309832741246e7; + + // levels energy + const double ha0 = 2.1471788680034824e10; + const double ha1 = 2.558764384495815e9; + const double ha2 = 5.323020344462938e8; + const double hb2 = 7.85178251911697e7; + + // Larmor frequency + double WL=0; + + + + complex E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies + + // density matrix elements which calculated via Hermitian property r_ij=conj(r_ji) + complex + r1301, + r1402, + r0903, + r1503, + r1004, + r1604, + r1105, + r0206, + r1406, + r0307, + r0907, + r1507, + r0408, + r1008, + r1608, + r1509, + r1610; + + + // inner use variables + double probability_v; // will be used as p(v) in Maxwell distribution + + ]]> + </globals> + <validation kind="run-time"/> <!--allows to put ranges as variables--> + <benchmark /> + <arguments> + <!-- Rabi frequency divided by 2 in [1/s] --> + <argument name="E1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" /> + <argument name="E2o" type="real" default_value="0.05*(2*M_PI*1e6)" /> + <argument name="E3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" /> + <argument name="E4o" type="real" default_value=".01*(2*M_PI*1e6)" /> + <!-- Fields detuning in [1/s] --> + <argument name="delta1" type="real" default_value="0.0" /> + <argument name="delta2" type="real" default_value="0.0" /> + <argument name="delta3" type="real" default_value="0.0" /> + <!--Pulse duration/width [s] --> + <argument name="Pwidth" type="real" default_value="0.1e-6" /> + <!-- Atom and cell properties --> + <!--Cell length [m] --> + <argument name="Lcell" type="real" default_value="1.5e-2" /> + <!--Density of atoms [1/m^3] --> + <argument name="Ndens" type="real" default_value="1e15" /> + <!--Atoms temperature [K] --> + <!--TODO: looks like Temperature > 10 K knocks solver, + I am guessing detunings are too large and thus it became a stiff equation--> + <!--! make sure it is not equal to zero!--> + <argument name="Temperature" type="real" default_value="5" /> + <!-- This will be executed after arguments/parameters are parsed --> + <!-- Read the code Luke: took me a while of reading the xmds2 sources to find it --> + <![CDATA[ + // Average sqrt(v^2) in Maxwell distribution for one dimension + if (Temperature == 0) + _LOG(_ERROR_LOG_LEVEL, "ERROR: Temperature should be >0 to provide range for Maxwell velocity distribution\n"); + v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); + // Maxwell distribution velocities range to take in account in [m/s] + // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0) + V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; + + // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)] + eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI; + // !FIXME over simplification: we should use relevant levels linewidths + eta1 = eta; + eta2 = eta; + eta3 = eta; + ]]> + </arguments> + <bing /> + <diagnostics /> + <fftw plan="estimate" threads="1" /> + <!-- I don't see any speed up on 6 core CPU even if use threads="6" --> + <!--<openmp />--> + <auto_vectorise /> + <halt_non_finite /> + </features> + + <!-- 'z', 't', and 'v' to have dimensions [m], [s], and [m/s] --> + <geometry> + <propagation_dimension> z </propagation_dimension> + <transverse_dimensions> + <!-- IMPORTANT: looks like having a lot of points in time helps with convergence. + I suspect that time spacing should be small enough to catch + all pulse harmonics and more importantly 1/dt should be larger than + the largest detuning (including Doppler shifts). + Unfortunately calculation time is proportional to lattice size + so we cannot just blindly increase it. + Some rules of thumb: + * lattice="1000" domain="(-1e-6, 1e-6)" + was good enough detunings up to 155 MHz (980 rad/s) notice that 1/dt=500 MHz + * lattice="10000" domain="(-1e-6, 1e-6)" + works for Doppler averaging in up to 400K for Rb when lasers are zero detuned + --> + <dimension name="t" lattice="10000" domain="(-1e-6, 1e-6)" /> + <dimension name="v" lattice="2" domain="(V_maxwell_min, V_maxwell_max)" /> + </transverse_dimensions> + </geometry> + + <!-- Rabi frequency --> + <vector name="E_field" type="complex" initial_space="t"> + <components>E1 E2 E3 E4</components> + <initialisation> + <![CDATA[ + // Initial (at starting 'z' position) electromagnetic field does not depend on detuning + // as well as time + E1=E1o; + E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) ); + E3=E3o; + E4=E4o; + ]]> + </initialisation> + </vector> + + <!--Maxwell distribution probability p(v)--> + <computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real"> + <components>probability_v</components> + <evaluation> + <![CDATA[ + // TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during the script lifetime since 'v' is fixed. + probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); + ]]> + </evaluation> + </computed_vector> + + <!--Maxwell distribution norm sum(p(v)) + Needed since we sum over the grid instead of true integral, + we also have finite cut off velocities--> + <computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real"> + <components>probability_v_norm</components> + <evaluation> + <dependencies basis="v">Maxwell_distribution_probabilities</dependencies> + <![CDATA[ + // TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during the script lifetime since 'v' is fixed. + probability_v_norm=probability_v; + ]]> + </evaluation> + </computed_vector> + + + <!-- Averaged across Maxwell distribution fields amplitudes --> + <computed_vector name="E_field_avgd" dimensions="t" type="complex"> + <components>E1a E2a E3a E4a</components> + <evaluation> + <dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies> + <![CDATA[ + double prob_v_normalized=probability_v/probability_v_norm; + E1a=E1*prob_v_normalized; + E2a=E2*prob_v_normalized; + E3a=E3*prob_v_normalized; + E4a=E4*prob_v_normalized; + ]]> + </evaluation> + </computed_vector> + + <!-- Averaged across Maxwell distribution density matrix components --> + <computed_vector name="density_matrix_averaged" dimensions="t" type="complex"> + <components> + r0101a + r0113a + r0202a + r0214a + r0303a + r0309a + r0315a + r0404a + r0410a + r0416a + r0505a + r0511a + r0602a + r0606a + r0614a + r0703a + r0707a + r0709a + r0715a + r0804a + r0808a + r0810a + r0816a + r0909a + r0915a + r1010a + r1016a + r1111a + r1313a + r1414a + r1515a + r1616a + </components> + <evaluation> + <dependencies basis="v">density_matrix Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies> + <![CDATA[ + double prob_v_normalized=probability_v/probability_v_norm; + + r0101a = r0101*prob_v_normalized; + r0113a = r0113*prob_v_normalized; + r0202a = r0202*prob_v_normalized; + r0214a = r0214*prob_v_normalized; + r0303a = r0303*prob_v_normalized; + r0309a = r0309*prob_v_normalized; + r0315a = r0315*prob_v_normalized; + r0404a = r0404*prob_v_normalized; + r0410a = r0410*prob_v_normalized; + r0416a = r0416*prob_v_normalized; + r0505a = r0505*prob_v_normalized; + r0511a = r0511*prob_v_normalized; + r0602a = r0602*prob_v_normalized; + r0606a = r0606*prob_v_normalized; + r0614a = r0614*prob_v_normalized; + r0703a = r0703*prob_v_normalized; + r0707a = r0707*prob_v_normalized; + r0709a = r0709*prob_v_normalized; + r0715a = r0715*prob_v_normalized; + r0804a = r0804*prob_v_normalized; + r0808a = r0808*prob_v_normalized; + r0810a = r0810*prob_v_normalized; + r0816a = r0816*prob_v_normalized; + r0909a = r0909*prob_v_normalized; + r0915a = r0915*prob_v_normalized; + r1010a = r1010*prob_v_normalized; + r1016a = r1016*prob_v_normalized; + r1111a = r1111*prob_v_normalized; + r1313a = r1313*prob_v_normalized; + r1414a = r1414*prob_v_normalized; + r1515a = r1515*prob_v_normalized; + r1616a = r1616*prob_v_normalized; + ]]> + </evaluation> + </computed_vector> + + + <vector name="density_matrix" type="complex" initial_space="t"> + <components> + r0101 + r0113 + r0202 + r0214 + r0303 + r0309 + r0315 + r0404 + r0410 + r0416 + r0505 + r0511 + r0602 + r0606 + r0614 + r0703 + r0707 + r0709 + r0715 + r0804 + r0808 + r0810 + r0816 + r0909 + r0915 + r1010 + r1016 + r1111 + r1313 + r1414 + r1515 + r1616 + </components> + <initialisation> + <!--This sets boundary condition at all times and left border of z (i.e. z=0)--> + <![CDATA[ + // Note: + // convergence is really slow if all populations concentrated at the bottom level |1> + // this is because if r11=1, everything else is 0 and then every small increment + // seems to be huge and adaptive solver makes smaller and smaller steps. + // As quick and dirty fix I reshuffle initial population + // so some of the population sits at the second ground level |2> + // TODO: Fix above. Make the equation of motion for r11 + // and express other level, let's say r44 + // through population normalization + r0101 = 0.125; + r0113 = 0; + r0202 = 0.125; + r0214 = 0; + r0303 = 0.125; + r0309 = 0; + r0315 = 0; + r0404 = 0.125; + r0410 = 0; + r0416 = 0; + r0505 = 0.125; + r0511 = 0; + r0602 = 0; + r0606 = 0.125; + r0614 = 0; + r0703 = 0; + r0707 = 0.125; + r0709 = 0; + r0715 = 0; + r0804 = 0; + r0808 = 0.125; + r0810 = 0; + r0816 = 0; + r0909 = 0; + r0915 = 0; + r1010 = 0; + r1016 = 0; + r1111 = 0; + r1313 = 0; + r1414 = 0; + r1515 = 0; + r1616 = 0; + ]]> + </initialisation> + </vector> + + <sequence> + <!--For this set of conditions ARK45 is faster than ARK89--> + <!--ARK45 is good for small detuning when all frequency like term are close to zero--> + <integrate algorithm="ARK45" tolerance="1e-5" interval="Lcell"> + <!--<integrate algorithm="SI" steps="200" interval="Lcell"> --> + <!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small--> + <!--When ARK45 works it about 3 times faster then RK4 with 1000 steps--> + <!--<integrate algorithm="RK4" steps="100" interval="1.5e-2">--> + <!--SIC algorithm seems to be much slower and needs fine 'z' step tuning and much finer time grid--> + <!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04--> + + <!--<integrate algorithm="SIC" interval="4e-2" steps="200">--> + + <samples>100</samples> + <!--<samples>100 100</samples>--> + <!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4--> + <!--<samples>100 100 100 100</samples>--> + <operators> + <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t"> + <integration_vectors>density_matrix</integration_vectors> + <dependencies>E_field_avgd</dependencies> + <boundary_condition kind="left"> + <!--This set boundary condition at all 'z' and left border of 't' (i.e. min(t))--> + <!-- + <![CDATA[ + r11 = 0; r22 = 1; r33 = 0; r44 = 0; + r12 = 0; r13 = 0; r14 = 0; + r23 = 0; r24 = 0; + r34 = 0; + printf("z= %g, t= %g\n", z, t); + ]]> + --> + </boundary_condition> + <![CDATA[ + E1ac = conj(E1a); + E2ac = conj(E2a); + E3ac = conj(E3a); + E4ac = conj(E4a); + + // Density matrix is Hermitian so we use r_ij=conj(r_ji) + + r1301 = conj(r0113); + r1402 = conj(r0214); + r0903 = conj(r0309); + r1503 = conj(r0315); + r1004 = conj(r0410); + r1604 = conj(r0416); + r1105 = conj(r0511); + r0206 = conj(r0602); + r1406 = conj(r0614); + r0307 = conj(r0703); + r0907 = conj(r0709); + r1507 = conj(r0715); + r0408 = conj(r0804); + r1008 = conj(r0810); + r1608 = conj(r0816); + r1509 = conj(r0915); + r1610 = conj(r1016); + + // Equations of motions according to Simon's mathematica code + dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (r1301*E4ac)/(4.*rt6)); + dr0113_dt = (-(gt*r0113) - (gt + g2)*r0113)/2. - i*(WL*r0113 - ((2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0113 + (r0101*E4ac)/(4.*rt6) - (r1313*E4ac)/(4.*rt6)); + dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (r1402*E4ac)/8.); + dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (r0206*E3ac)/(8.*rt3) + (r0202*E4ac)/8. - (r1414*E4ac)/8.); + dr0303_dt = gt/8. - gt*r0303 + (g1*r0909)/12. + (g1*r1010)/3. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1515)/4. - i*((r0309*E1a)/(4.*rt6) + (r0315*E4a)/8. - (r0903*E1ac)/(4.*rt6) - (r1503*E4ac)/8.); + dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*E1ac)/(4.*rt6) - (r0909*E1ac)/(4.*rt6) - (r0307*E2ac)/(4.*rt6) - (r1509*E4ac)/8.); + dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*E1ac)/(4.*rt6) - (r0307*E3ac)/8. + (r0303*E4ac)/8. - (r1515*E4ac)/8.); + dr0404_dt = gt/8. - gt*r0404 + (g1*r1010)/4. + (g1*r1111)/4. + (g2*r1414)/4. + (g2*r1515)/12. + (g2*r1616)/6. - i*((r0410*E1a)/(4.*rt2) + (r0416*E4a)/(4.*rt6) - (r1004*E1ac)/(4.*rt2) - (r1604*E4ac)/(4.*rt6)); + dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*E1ac)/(4.*rt2) - (r1010*E1ac)/(4.*rt2) - (r0408*E2ac)/(4.*rt6) - (r1610*E4ac)/(4.*rt6)); + dr0416_dt = (-(gt*r0416) - (gt + g2)*r0416)/2. - i*(-(WL*r0416)/2. - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0416 - (r1016*E1ac)/(4.*rt2) - (r0408*E3ac)/(4.*rt2) + (r0404*E4ac)/(4.*rt6) - (r1616*E4ac)/(4.*rt6)); + dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (r1105*E1ac)/4.); + dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*E1ac)/4. - (r1111*E1ac)/4.); + dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (r1402*E3ac)/(8.*rt3)); + dr0606_dt = gt/8. - gt*r0606 + (g1*r0909)/12. + (g1*r1010)/12. + (g2*r1313)/4. + (g2*r1414)/12. - i*(-(r0614*E3a)/(8.*rt3) + (r1406*E3ac)/(8.*rt3)); + dr0614_dt = (-(gt*r0614) - (gt + g2)*r0614)/2. - i*((-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0614 - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0614 - (r0606*E3ac)/(8.*rt3) + (r1414*E3ac)/(8.*rt3) + (r0602*E4ac)/8.); + dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (r0903*E2ac)/(4.*rt6) + (r1503*E3ac)/8.); + dr0707_dt = gt/8. - gt*r0707 + (g1*r0909)/12. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1414)/3. + (g2*r1515)/4. - i*(-(r0709*E2a)/(4.*rt6) - (r0715*E3a)/8. + (r0907*E2ac)/(4.*rt6) + (r1507*E3ac)/8.); + dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*E1ac)/(4.*rt6) - (r0707*E2ac)/(4.*rt6) + (r0909*E2ac)/(4.*rt6) + (r1509*E3ac)/8.); + dr0715_dt = (-(gt*r0715) - (gt + g2)*r0715)/2. - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0715 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0715 + (r0915*E2ac)/(4.*rt6) - (r0707*E3ac)/8. + (r1515*E3ac)/8. + (r0703*E4ac)/8.); + dr0804_dt = -(gt*r0804) - i*((WL*r0804)/2. + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0804 + (r0810*E1a)/(4.*rt2) + (r0816*E4a)/(4.*rt6) + (r1004*E2ac)/(4.*rt6) + (r1604*E3ac)/(4.*rt2)); + dr0808_dt = gt/8. - gt*r0808 + (g1*r1010)/12. + (g1*r1111)/12. + (g2*r1414)/12. + (g2*r1515)/4. + (g2*r1616)/2. - i*(-(r0810*E2a)/(4.*rt6) - (r0816*E3a)/(4.*rt2) + (r1008*E2ac)/(4.*rt6) + (r1608*E3ac)/(4.*rt2)); + dr0810_dt = (-(gt*r0810) - (gt + g1)*r0810)/2. - i*((delta1 + v*Kvec1)*r0810 + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0810 + (r0804*E1ac)/(4.*rt2) - (r0808*E2ac)/(4.*rt6) + (r1010*E2ac)/(4.*rt6) + (r1610*E3ac)/(4.*rt2)); + dr0816_dt = (-(gt*r0816) - (gt + g2)*r0816)/2. - i*((WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0816 - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0816 + (r1016*E2ac)/(4.*rt6) - (r0808*E3ac)/(4.*rt2) + (r1616*E3ac)/(4.*rt2) + (r0804*E4ac)/(4.*rt6)); + dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (r0903*E1ac)/(4.*rt6) - (r0907*E2ac)/(4.*rt6)); + dr0915_dt = (-((gt + g1)*r0915) - (gt + g2)*r0915)/2. - i*((-WL/6. - delta1 - v*Kvec1)*r0915 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0915 - (r0315*E1a)/(4.*rt6) + (r0715*E2a)/(4.*rt6) - (r0907*E3ac)/8. + (r0903*E4ac)/8.); + dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (r1004*E1ac)/(4.*rt2) - (r1008*E2ac)/(4.*rt6)); + dr1016_dt = (-((gt + g1)*r1016) - (gt + g2)*r1016)/2. - i*(-((delta1 + v*Kvec1)*r1016) - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1016 - (r0416*E1a)/(4.*rt2) + (r0816*E2a)/(4.*rt6) - (r1008*E3ac)/(4.*rt2) + (r1004*E4ac)/(4.*rt6)); + dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (r1105*E1ac)/4.); + dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (r1301*E4ac)/(4.*rt6)); + dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (r1406*E3ac)/(8.*rt3) + (r1402*E4ac)/8.); + dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (r1507*E3ac)/8. + (r1503*E4ac)/8.); + dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (r1608*E3ac)/(4.*rt2) + (r1604*E4ac)/(4.*rt6)); + ]]> + </operator> + <!-- + According to xmds2 docs operator kind="ip" should be faster + but our codes runs about 5% to 10% slower with it. + Maybe because we very close to the stiff condition so I use "ex" kind + <operator kind="ip" constant="yes"> + --> + <operator kind="ex" constant="yes" type="imaginary"> + <operator_names>Lt</operator_names> + <![CDATA[ + Lt = -i/c*kt; + ]]> + </operator> + <integration_vectors>E_field</integration_vectors> + <dependencies>density_matrix</dependencies> + <![CDATA[ + // Density matrix is Hermitian so we use r_ij=conj(r_ji) + r1301 = conj(r0113); + r1402 = conj(r0214); + r0903 = conj(r0309); + r1503 = conj(r0315); + r1004 = conj(r0410); + r1604 = conj(r0416); + r1105 = conj(r0511); + r0206 = conj(r0602); + r1406 = conj(r0614); + r0307 = conj(r0703); + r0907 = conj(r0709); + r1507 = conj(r0715); + r0408 = conj(r0804); + r1008 = conj(r0810); + r1608 = conj(r0816); + r1509 = conj(r0915); + r1610 = conj(r1016); + + dE1_dz = 0.16666666666666666*i*eta1*(2.449489742783178*r0903 + 4.242640687119286*r1004 + 6.*r1105) - Lt[E1]; + dE2_dz = -0.4082482904638631*i*eta1*(r0907 + r1008) - Lt[E2]; + dE3_dz = -0.3333333333333333*i*eta2*(1.7320508075688772*r1406 + 3.*r1507 + 4.242640687119286*r1608) - Lt[E3]; + dE4_dz = (i*eta2*(2.449489742783178*r1301 + 3*r1402 + 3*r1503 + 2.449489742783178*r1604))/3. - Lt[E4]; + ]]> + </operators> + </integrate> + </sequence> + + + + <!-- The output to generate --> + <output format="binary" filename="realistic_Rb.xsil"> + <group> + <sampling basis="t(1000) " initial_sample="yes"> + <dependencies>E_field_avgd</dependencies> + <moments>I1_out I2_out I3_out I4_out</moments> + <![CDATA[ + I1_out = mod2(E1a); + I2_out = mod2(E2a); + I3_out = mod2(E3a); + I4_out = mod2(E4a); + ]]> + </sampling> + </group> + + <!-- + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>density_matrix_averaged</dependencies> + <moments> + r11_out r22_out r33_out r44_out + r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out + r23_re_out r23_im_out r24_re_out r24_im_out + r34_re_out r34_im_out + </moments> + <![CDATA[ + // populations output + r11_out = r11a.Re(); + r22_out = r22a.Re(); + r33_out = r33a.Re(); + r44_out = r44a.Re(); + // coherences output + r12_re_out = r12a.Re(); + r12_im_out = r12a.Im(); + r13_re_out = r13a.Re(); + r13_im_out = r13a.Im(); + r14_re_out = r14a.Re(); + r14_im_out = r14a.Im(); + r23_re_out = r23a.Re(); + r23_im_out = r23a.Im(); + r24_re_out = r24a.Re(); + r24_im_out = r24a.Im(); + r34_re_out = r34a.Re(); + r34_im_out = r34a.Im(); + ]]> + </sampling> + </group> + --> + + <!-- use the following two groups only for debuging + otherwise they are quite useless and have to much information + in 3D space (z,t,v) --> + <!-- + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>E_field</dependencies> + <moments>I1_out_v I2_out_v I3_out_v I4_out_v</moments> + <![CDATA[ + // light field intensity distribution in velocity subgroups + I1_out_v = mod2(E1); + I2_out_v = mod2(E2); + I3_out_v = mod2(E3); + I4_out_v = mod2(E4); + ]]> + </sampling> + </group> + + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>density_matrix</dependencies> + <moments> + r11_out_v r22_out_v r33_out_v r44_out_v + r12_re_out_v r12_im_out_v r13_re_out_v r13_im_out_v r14_re_out_v r14_im_out_v + r23_re_out_v r23_im_out_v r24_re_out_v r24_im_out_v + r34_re_out_v r34_im_out_v + </moments> + <![CDATA[ + // density matrix distribution in velocity subgroups + // populations output + r11_out_v = r11.Re(); + r22_out_v = r22.Re(); + r33_out_v = r33.Re(); + r44_out_v = r44.Re(); + // coherences output + r12_re_out_v = r12.Re(); + r12_im_out_v = r12.Im(); + r13_re_out_v = r13.Re(); + r13_im_out_v = r13.Im(); + r14_re_out_v = r14.Re(); + r14_im_out_v = r14.Im(); + r23_re_out_v = r23.Re(); + r23_im_out_v = r23.Im(); + r24_re_out_v = r24.Re(); + r24_im_out_v = r24.Im(); + r34_re_out_v = r34.Re(); + r34_im_out_v = r34.Im(); + ]]> + </sampling> + </group> + --> + + </output> + +</simulation> + +<!-- +vim: ts=2 sw=2 foldmethod=indent: +--> diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil new file mode 100644 index 0000000..409374d --- /dev/null +++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected.xsil @@ -0,0 +1,699 @@ +<?xml version="1.0"?> +<simulation xmds-version="2"> + <testing> + <arguments> --Ndens=1e15 --Lcell=10.0e-2 --Temperature=1e-3 --Pwidth=0.4e-6 --delta1=0 --delta2=0 --delta3=0 --E1o=0 --E2o=1e2 --E3o=0 --E4o=0</arguments> + <xsil_file name="realistic_Rb.xsil" expected="realistic_Rb_expected.xsil" absolute_tolerance="1e-7" relative_tolerance="1e-5"> + <moment_group number="0" absolute_tolerance="1e-7" relative_tolerance="1e-6" /> + </xsil_file> + </testing> + + <name>realistic_Rb</name> + + <author>Eugeniy Mikhailov</author> + <description> + License GPL. + + Solving simplified Rb atom model + with fields propagation along spatial axis Z + with Doppler broadening. + + + We assume four-wave mixing condition when w3-w4=w2-w1 i.e. fields E3 and E4 drive the same + resonance as fields E2 and E1. + + + * --------------- | F=1, 2P_3/2 > + * \ \ + * \ E3_r \ -------- | F=2, 2P_+1/2 > + * \ E4_r \ / \ + * \ \ / E2_l \ + * \ / \ E1_l + * | F=2, 2S_1/2 > -------------- \ + * \ \ + * \ \ + * ------------- | F=1, 2S_1/2 > + * + + + We are solving + dE/dz+(1/c)*dE/dt=i*eta*rho_ij, where j level is higher then i. + Note that E is actually a Rabi frequency of electromagnetic field not the EM field + in xmds terms it looks like + dE_dz = i*eta*rhoij - 1/c*L[E], here we moved t dependence to Fourier space + + VERY IMPORTANT: all Rabi frequency should be given in [1/s], if you want to + normalize it to something else look drho/dt equation. + No need to renormalizes eta as long as its express through + the upper level decay rate in the same units as Rabi frequency. + </description> + + <features> + <globals> + <![CDATA[ + // Some numerical constants + const double pi = M_PI; + // proportional to splitting ratios sqrt(6) , sqrt(3), sqrt(2) + const double rt6 = 2.449489742783178; + const double rt3 = 1.7320508075688772; + const double rt2 = 1.4142135623730951; + + + const double c=3.e8; + const double k_boltzmann= 1.3806505e-23; // Boltzmann knostant in [J/K] + const double lambda=794.7e-9; //wavelength in m + // Fields k-vector + const double Kvec = 2*M_PI/lambda; + // Simplified k-vectors + const double Kvec1 = Kvec, Kvec2=Kvec, Kvec3=Kvec; + + const double Gamma_super=6*(2*M_PI*1e6); // characteristic decay rate of upper level used for eta calculations expressed in [1/s] + // eta will be calculated in the <arguments> section + double eta = 0; // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)] + double eta1=0, eta2=0, eta3=0; + + // --------- Atom and cell properties ------------------------- + // range of Maxwell distribution atomic velocities + const double mass = (86.909180527 * 1.660538921e-27); // atom mass in [kg] + // above mass expression is written as (expression is isotopic_mass * atomic_mass_unit) + + // Average sqrt(v^2) in Maxwell distribution for one dimension + // Maxwell related parameters will be calculated in <arguments> section + double v_thermal_averaged=0; + // Maxwell distribution velocities range to take in account in [m/s] + double V_maxwell_min = 0, V_maxwell_max = 0; + + // repopulation rate (atoms flying in/out the laser beam) in [1/s] + const double gt=0.01 *(2*M_PI*1e6); + + // Natural linewidth of j's level in [1/s] + const double g1 = 3.612847284945266e7; + const double g2 = 3.8117309832741246e7; + + // levels energy + const double ha0 = 2.1471788680034824e10; + const double ha1 = 2.558764384495815e9; + const double ha2 = 5.323020344462938e8; + const double hb2 = 7.85178251911697e7; + + // Larmor frequency + double WL=0; + + + + complex E1ac, E2ac, E3ac, E4ac; // Complex conjugated Rabi frequencies + + // density matrix elements which calculated via Hermitian property r_ij=conj(r_ji) + complex + r1301, + r1402, + r0903, + r1503, + r1004, + r1604, + r1105, + r0206, + r1406, + r0307, + r0907, + r1507, + r0408, + r1008, + r1608, + r1509, + r1610; + + + // inner use variables + double probability_v; // will be used as p(v) in Maxwell distribution + + ]]> + </globals> + <validation kind="run-time"/> <!--allows to put ranges as variables--> + <benchmark /> + <arguments> + <!-- Rabi frequency divided by 2 in [1/s] --> + <argument name="E1o" type="real" default_value="2*1.5*(2*M_PI*1e6)" /> + <argument name="E2o" type="real" default_value="0.05*(2*M_PI*1e6)" /> + <argument name="E3o" type="real" default_value="2*3.0*(2*M_PI*1e6)" /> + <argument name="E4o" type="real" default_value=".01*(2*M_PI*1e6)" /> + <!-- Fields detuning in [1/s] --> + <argument name="delta1" type="real" default_value="0.0" /> + <argument name="delta2" type="real" default_value="0.0" /> + <argument name="delta3" type="real" default_value="0.0" /> + <!--Pulse duration/width [s] --> + <argument name="Pwidth" type="real" default_value="0.1e-6" /> + <!-- Atom and cell properties --> + <!--Cell length [m] --> + <argument name="Lcell" type="real" default_value="1.5e-2" /> + <!--Density of atoms [1/m^3] --> + <argument name="Ndens" type="real" default_value="1e15" /> + <!--Atoms temperature [K] --> + <!--TODO: looks like Temperature > 10 K knocks solver, + I am guessing detunings are too large and thus it became a stiff equation--> + <!--! make sure it is not equal to zero!--> + <argument name="Temperature" type="real" default_value="5" /> + <!-- This will be executed after arguments/parameters are parsed --> + <!-- Read the code Luke: took me a while of reading the xmds2 sources to find it --> + <![CDATA[ + // Average sqrt(v^2) in Maxwell distribution for one dimension + if (Temperature == 0) + _LOG(_ERROR_LOG_LEVEL, "ERROR: Temperature should be >0 to provide range for Maxwell velocity distribution\n"); + v_thermal_averaged=sqrt(k_boltzmann*Temperature/mass); + // Maxwell distribution velocities range to take in account in [m/s] + // there is almost zero probability for higher velocity p(4*v_av) = 3.3e-04 * p(0) + V_maxwell_min = -4*v_thermal_averaged; V_maxwell_max = -V_maxwell_min; + + // eta constant in the wave equation for Rabi frequency. Units are [1/(m s)] + eta = 3*lambda*lambda*Ndens*Gamma_super/8.0/M_PI; + // !FIXME over simplification: we should use relevant levels linewidths + eta1 = eta; + eta2 = eta; + eta3 = eta; + ]]> + </arguments> + <bing /> + <diagnostics /> + <fftw plan="estimate" threads="1" /> + <!-- I don't see any speed up on 6 core CPU even if use threads="6" --> + <!--<openmp />--> + <auto_vectorise /> + <halt_non_finite /> + </features> + + <!-- 'z', 't', and 'v' to have dimensions [m], [s], and [m/s] --> + <geometry> + <propagation_dimension> z </propagation_dimension> + <transverse_dimensions> + <!-- IMPORTANT: looks like having a lot of points in time helps with convergence. + I suspect that time spacing should be small enough to catch + all pulse harmonics and more importantly 1/dt should be larger than + the largest detuning (including Doppler shifts). + Unfortunately calculation time is proportional to lattice size + so we cannot just blindly increase it. + Some rules of thumb: + * lattice="1000" domain="(-1e-6, 1e-6)" + was good enough detunings up to 155 MHz (980 rad/s) notice that 1/dt=500 MHz + * lattice="10000" domain="(-1e-6, 1e-6)" + works for Doppler averaging in up to 400K for Rb when lasers are zero detuned + --> + <dimension name="t" lattice="10000" domain="(-1e-6, 1e-6)" /> + <dimension name="v" lattice="2" domain="(V_maxwell_min, V_maxwell_max)" /> + </transverse_dimensions> + </geometry> + + <!-- Rabi frequency --> + <vector name="E_field" type="complex" initial_space="t"> + <components>E1 E2 E3 E4</components> + <initialisation> + <![CDATA[ + // Initial (at starting 'z' position) electromagnetic field does not depend on detuning + // as well as time + E1=E1o; + E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) ); + E3=E3o; + E4=E4o; + ]]> + </initialisation> + </vector> + + <!--Maxwell distribution probability p(v)--> + <computed_vector name="Maxwell_distribution_probabilities" dimensions="v" type="real"> + <components>probability_v</components> + <evaluation> + <![CDATA[ + // TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during the script lifetime since 'v' is fixed. + probability_v=1.0/(v_thermal_averaged*sqrt(2*M_PI)) * exp( - mod2(v/v_thermal_averaged)/2.0 ); + ]]> + </evaluation> + </computed_vector> + + <!--Maxwell distribution norm sum(p(v)) + Needed since we sum over the grid instead of true integral, + we also have finite cut off velocities--> + <computed_vector name="Maxwell_distribution_probabilities_norm" dimensions="" type="real"> + <components>probability_v_norm</components> + <evaluation> + <dependencies basis="v">Maxwell_distribution_probabilities</dependencies> + <![CDATA[ + // TODO: move to the global space/function. This reevaluated many times since it called from dependency requests but it never changes during the script lifetime since 'v' is fixed. + probability_v_norm=probability_v; + ]]> + </evaluation> + </computed_vector> + + + <!-- Averaged across Maxwell distribution fields amplitudes --> + <computed_vector name="E_field_avgd" dimensions="t" type="complex"> + <components>E1a E2a E3a E4a</components> + <evaluation> + <dependencies basis="v">E_field Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies> + <![CDATA[ + double prob_v_normalized=probability_v/probability_v_norm; + E1a=E1*prob_v_normalized; + E2a=E2*prob_v_normalized; + E3a=E3*prob_v_normalized; + E4a=E4*prob_v_normalized; + ]]> + </evaluation> + </computed_vector> + + <!-- Averaged across Maxwell distribution density matrix components --> + <computed_vector name="density_matrix_averaged" dimensions="t" type="complex"> + <components> + r0101a + r0113a + r0202a + r0214a + r0303a + r0309a + r0315a + r0404a + r0410a + r0416a + r0505a + r0511a + r0602a + r0606a + r0614a + r0703a + r0707a + r0709a + r0715a + r0804a + r0808a + r0810a + r0816a + r0909a + r0915a + r1010a + r1016a + r1111a + r1313a + r1414a + r1515a + r1616a + </components> + <evaluation> + <dependencies basis="v">density_matrix Maxwell_distribution_probabilities Maxwell_distribution_probabilities_norm</dependencies> + <![CDATA[ + double prob_v_normalized=probability_v/probability_v_norm; + + r0101a = r0101*prob_v_normalized; + r0113a = r0113*prob_v_normalized; + r0202a = r0202*prob_v_normalized; + r0214a = r0214*prob_v_normalized; + r0303a = r0303*prob_v_normalized; + r0309a = r0309*prob_v_normalized; + r0315a = r0315*prob_v_normalized; + r0404a = r0404*prob_v_normalized; + r0410a = r0410*prob_v_normalized; + r0416a = r0416*prob_v_normalized; + r0505a = r0505*prob_v_normalized; + r0511a = r0511*prob_v_normalized; + r0602a = r0602*prob_v_normalized; + r0606a = r0606*prob_v_normalized; + r0614a = r0614*prob_v_normalized; + r0703a = r0703*prob_v_normalized; + r0707a = r0707*prob_v_normalized; + r0709a = r0709*prob_v_normalized; + r0715a = r0715*prob_v_normalized; + r0804a = r0804*prob_v_normalized; + r0808a = r0808*prob_v_normalized; + r0810a = r0810*prob_v_normalized; + r0816a = r0816*prob_v_normalized; + r0909a = r0909*prob_v_normalized; + r0915a = r0915*prob_v_normalized; + r1010a = r1010*prob_v_normalized; + r1016a = r1016*prob_v_normalized; + r1111a = r1111*prob_v_normalized; + r1313a = r1313*prob_v_normalized; + r1414a = r1414*prob_v_normalized; + r1515a = r1515*prob_v_normalized; + r1616a = r1616*prob_v_normalized; + ]]> + </evaluation> + </computed_vector> + + + <vector name="density_matrix" type="complex" initial_space="t"> + <components> + r0101 + r0113 + r0202 + r0214 + r0303 + r0309 + r0315 + r0404 + r0410 + r0416 + r0505 + r0511 + r0602 + r0606 + r0614 + r0703 + r0707 + r0709 + r0715 + r0804 + r0808 + r0810 + r0816 + r0909 + r0915 + r1010 + r1016 + r1111 + r1313 + r1414 + r1515 + r1616 + </components> + <initialisation> + <!--This sets boundary condition at all times and left border of z (i.e. z=0)--> + <![CDATA[ + // Note: + // convergence is really slow if all populations concentrated at the bottom level |1> + // this is because if r11=1, everything else is 0 and then every small increment + // seems to be huge and adaptive solver makes smaller and smaller steps. + // As quick and dirty fix I reshuffle initial population + // so some of the population sits at the second ground level |2> + // TODO: Fix above. Make the equation of motion for r11 + // and express other level, let's say r44 + // through population normalization + r0101 = 0.125; + r0113 = 0; + r0202 = 0.125; + r0214 = 0; + r0303 = 0.125; + r0309 = 0; + r0315 = 0; + r0404 = 0.125; + r0410 = 0; + r0416 = 0; + r0505 = 0.125; + r0511 = 0; + r0602 = 0; + r0606 = 0.125; + r0614 = 0; + r0703 = 0; + r0707 = 0.125; + r0709 = 0; + r0715 = 0; + r0804 = 0; + r0808 = 0.125; + r0810 = 0; + r0816 = 0; + r0909 = 0; + r0915 = 0; + r1010 = 0; + r1016 = 0; + r1111 = 0; + r1313 = 0; + r1414 = 0; + r1515 = 0; + r1616 = 0; + ]]> + </initialisation> + </vector> + + <sequence> + <!--For this set of conditions ARK45 is faster than ARK89--> + <!--ARK45 is good for small detuning when all frequency like term are close to zero--> + <integrate algorithm="ARK45" tolerance="1e-5" interval="Lcell"> + <!--<integrate algorithm="SI" steps="200" interval="Lcell"> --> + <!--RK4 is good for large detunings when frequency like term are big, it does not try to be too smart about adaptive step which ARK seems to make too small--> + <!--When ARK45 works it about 3 times faster then RK4 with 1000 steps--> + <!--<integrate algorithm="RK4" steps="100" interval="1.5e-2">--> + <!--SIC algorithm seems to be much slower and needs fine 'z' step tuning and much finer time grid--> + <!--For example I had to quadruple the time grid from 1000 to 4000 when increased z distance from 0.02 to 0.04--> + + <!--<integrate algorithm="SIC" interval="4e-2" steps="200">--> + + <samples>100</samples> + <!--<samples>100 100</samples>--> + <!--Use the next line for debuging to see velocity dependence. Uncomment/switch on output groups 3,4--> + <!--<samples>100 100 100 100</samples>--> + <operators> + <operator kind="cross_propagation" algorithm="SI" propagation_dimension="t"> + <integration_vectors>density_matrix</integration_vectors> + <dependencies>E_field_avgd</dependencies> + <boundary_condition kind="left"> + <!--This set boundary condition at all 'z' and left border of 't' (i.e. min(t))--> + <!-- + <![CDATA[ + r11 = 0; r22 = 1; r33 = 0; r44 = 0; + r12 = 0; r13 = 0; r14 = 0; + r23 = 0; r24 = 0; + r34 = 0; + printf("z= %g, t= %g\n", z, t); + ]]> + --> + </boundary_condition> + <![CDATA[ + E1ac = conj(E1a); + E2ac = conj(E2a); + E3ac = conj(E3a); + E4ac = conj(E4a); + + // Density matrix is Hermitian so we use r_ij=conj(r_ji) + + r1301 = conj(r0113); + r1402 = conj(r0214); + r0903 = conj(r0309); + r1503 = conj(r0315); + r1004 = conj(r0410); + r1604 = conj(r0416); + r1105 = conj(r0511); + r0206 = conj(r0602); + r1406 = conj(r0614); + r0307 = conj(r0703); + r0907 = conj(r0709); + r1507 = conj(r0715); + r0408 = conj(r0804); + r1008 = conj(r0810); + r1608 = conj(r0816); + r1509 = conj(r0915); + r1610 = conj(r1016); + + // Equations of motions according to Simon's mathematica code + dr0101_dt = gt/8. - gt*r0101 + (g1*r0909)/2. + (g2*r1313)/6. - i*((r0113*E4a)/(4.*rt6) - (r1301*E4ac)/(4.*rt6)); + dr0113_dt = (-(gt*r0113) - (gt + g2)*r0113)/2. - i*(WL*r0113 - ((2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0113 + (r0101*E4ac)/(4.*rt6) - (r1313*E4ac)/(4.*rt6)); + dr0202_dt = gt/8. - gt*r0202 + (g1*r0909)/4. + (g1*r1010)/4. + (g2*r1313)/12. + (g2*r1414)/4. - i*((r0214*E4a)/8. - (r1402*E4ac)/8.); + dr0214_dt = (-(gt*r0214) - (gt + g2)*r0214)/2. - i*((WL*r0214)/2. - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0214 - (r0206*E3ac)/(8.*rt3) + (r0202*E4ac)/8. - (r1414*E4ac)/8.); + dr0303_dt = gt/8. - gt*r0303 + (g1*r0909)/12. + (g1*r1010)/3. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1515)/4. - i*((r0309*E1a)/(4.*rt6) + (r0315*E4a)/8. - (r0903*E1ac)/(4.*rt6) - (r1503*E4ac)/8.); + dr0309_dt = (-(gt*r0309) - (gt + g1)*r0309)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0309) + (r0303*E1ac)/(4.*rt6) - (r0909*E1ac)/(4.*rt6) - (r0307*E2ac)/(4.*rt6) - (r1509*E4ac)/8.); + dr0315_dt = (-(gt*r0315) - (gt + g2)*r0315)/2. - i*(-(((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0315) - (r0915*E1ac)/(4.*rt6) - (r0307*E3ac)/8. + (r0303*E4ac)/8. - (r1515*E4ac)/8.); + dr0404_dt = gt/8. - gt*r0404 + (g1*r1010)/4. + (g1*r1111)/4. + (g2*r1414)/4. + (g2*r1515)/12. + (g2*r1616)/6. - i*((r0410*E1a)/(4.*rt2) + (r0416*E4a)/(4.*rt6) - (r1004*E1ac)/(4.*rt2) - (r1604*E4ac)/(4.*rt6)); + dr0410_dt = (-(gt*r0410) - (gt + g1)*r0410)/2. - i*(-(WL*r0410)/2. + (delta1 + v*Kvec1)*r0410 + (r0404*E1ac)/(4.*rt2) - (r1010*E1ac)/(4.*rt2) - (r0408*E2ac)/(4.*rt6) - (r1610*E4ac)/(4.*rt6)); + dr0416_dt = (-(gt*r0416) - (gt + g2)*r0416)/2. - i*(-(WL*r0416)/2. - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0416 - (r1016*E1ac)/(4.*rt2) - (r0408*E3ac)/(4.*rt2) + (r0404*E4ac)/(4.*rt6) - (r1616*E4ac)/(4.*rt6)); + dr0505_dt = gt/8. - gt*r0505 + (g1*r1111)/2. + (g2*r1515)/6. + (g2*r1616)/3. - i*((r0511*E1a)/4. - (r1105*E1ac)/4.); + dr0511_dt = (-(gt*r0511) - (gt + g1)*r0511)/2. - i*(-(WL*r0511) - (WL/6. - delta1 - v*Kvec1)*r0511 + (r0505*E1ac)/4. - (r1111*E1ac)/4.); + dr0602_dt = -(gt*r0602) - i*(-(WL*r0602)/2. + (-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0602 + (r0614*E4a)/8. + (r1402*E3ac)/(8.*rt3)); + dr0606_dt = gt/8. - gt*r0606 + (g1*r0909)/12. + (g1*r1010)/12. + (g2*r1313)/4. + (g2*r1414)/12. - i*(-(r0614*E3a)/(8.*rt3) + (r1406*E3ac)/(8.*rt3)); + dr0614_dt = (-(gt*r0614) - (gt + g2)*r0614)/2. - i*((-WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0614 - (-delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0614 - (r0606*E3ac)/(8.*rt3) + (r1414*E3ac)/(8.*rt3) + (r0602*E4ac)/8.); + dr0703_dt = -(gt*r0703) - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0703 + (r0709*E1a)/(4.*rt6) + (r0715*E4a)/8. + (r0903*E2ac)/(4.*rt6) + (r1503*E3ac)/8.); + dr0707_dt = gt/8. - gt*r0707 + (g1*r0909)/12. + (g1*r1111)/12. + (g2*r1313)/4. + (g2*r1414)/3. + (g2*r1515)/4. - i*(-(r0709*E2a)/(4.*rt6) - (r0715*E3a)/8. + (r0907*E2ac)/(4.*rt6) + (r1507*E3ac)/8.); + dr0709_dt = (-(gt*r0709) - (gt + g1)*r0709)/2. - i*(-((-WL/6. - delta1 - v*Kvec1)*r0709) + (-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0709 + (r0703*E1ac)/(4.*rt6) - (r0707*E2ac)/(4.*rt6) + (r0909*E2ac)/(4.*rt6) + (r1509*E3ac)/8.); + dr0715_dt = (-(gt*r0715) - (gt + g2)*r0715)/2. - i*((-delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0715 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0715 + (r0915*E2ac)/(4.*rt6) - (r0707*E3ac)/8. + (r1515*E3ac)/8. + (r0703*E4ac)/8.); + dr0804_dt = -(gt*r0804) - i*((WL*r0804)/2. + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0804 + (r0810*E1a)/(4.*rt2) + (r0816*E4a)/(4.*rt6) + (r1004*E2ac)/(4.*rt6) + (r1604*E3ac)/(4.*rt2)); + dr0808_dt = gt/8. - gt*r0808 + (g1*r1010)/12. + (g1*r1111)/12. + (g2*r1414)/12. + (g2*r1515)/4. + (g2*r1616)/2. - i*(-(r0810*E2a)/(4.*rt6) - (r0816*E3a)/(4.*rt2) + (r1008*E2ac)/(4.*rt6) + (r1608*E3ac)/(4.*rt2)); + dr0810_dt = (-(gt*r0810) - (gt + g1)*r0810)/2. - i*((delta1 + v*Kvec1)*r0810 + (WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0810 + (r0804*E1ac)/(4.*rt2) - (r0808*E2ac)/(4.*rt6) + (r1010*E2ac)/(4.*rt6) + (r1610*E3ac)/(4.*rt2)); + dr0816_dt = (-(gt*r0816) - (gt + g2)*r0816)/2. - i*((WL/2. - delta1 + delta2 - v*Kvec1 + v*Kvec2)*r0816 - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0816 + (r1016*E2ac)/(4.*rt6) - (r0808*E3ac)/(4.*rt2) + (r1616*E3ac)/(4.*rt2) + (r0804*E4ac)/(4.*rt6)); + dr0909_dt = -((gt + g1)*r0909) - i*(-(r0309*E1a)/(4.*rt6) + (r0709*E2a)/(4.*rt6) + (r0903*E1ac)/(4.*rt6) - (r0907*E2ac)/(4.*rt6)); + dr0915_dt = (-((gt + g1)*r0915) - (gt + g2)*r0915)/2. - i*((-WL/6. - delta1 - v*Kvec1)*r0915 - ((-2*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r0915 - (r0315*E1a)/(4.*rt6) + (r0715*E2a)/(4.*rt6) - (r0907*E3ac)/8. + (r0903*E4ac)/8.); + dr1010_dt = -((gt + g1)*r1010) - i*(-(r0410*E1a)/(4.*rt2) + (r0810*E2a)/(4.*rt6) + (r1004*E1ac)/(4.*rt2) - (r1008*E2ac)/(4.*rt6)); + dr1016_dt = (-((gt + g1)*r1016) - (gt + g2)*r1016)/2. - i*(-((delta1 + v*Kvec1)*r1016) - ((-4*WL)/3. - delta1 + delta2 - delta3 - v*Kvec1 + v*Kvec2 - v*Kvec3)*r1016 - (r0416*E1a)/(4.*rt2) + (r0816*E2a)/(4.*rt6) - (r1008*E3ac)/(4.*rt2) + (r1004*E4ac)/(4.*rt6)); + dr1111_dt = -((gt + g1)*r1111) - i*(-(r0511*E1a)/4. + (r1105*E1ac)/4.); + dr1313_dt = -((gt + g2)*r1313) - i*(-(r0113*E4a)/(4.*rt6) + (r1301*E4ac)/(4.*rt6)); + dr1414_dt = -((gt + g2)*r1414) - i*((r0614*E3a)/(8.*rt3) - (r0214*E4a)/8. - (r1406*E3ac)/(8.*rt3) + (r1402*E4ac)/8.); + dr1515_dt = -((gt + g2)*r1515) - i*((r0715*E3a)/8. - (r0315*E4a)/8. - (r1507*E3ac)/8. + (r1503*E4ac)/8.); + dr1616_dt = -((gt + g2)*r1616) - i*((r0816*E3a)/(4.*rt2) - (r0416*E4a)/(4.*rt6) - (r1608*E3ac)/(4.*rt2) + (r1604*E4ac)/(4.*rt6)); + ]]> + </operator> + <!-- + According to xmds2 docs operator kind="ip" should be faster + but our codes runs about 5% to 10% slower with it. + Maybe because we very close to the stiff condition so I use "ex" kind + <operator kind="ip" constant="yes"> + --> + <operator kind="ex" constant="yes" type="imaginary"> + <operator_names>Lt</operator_names> + <![CDATA[ + Lt = -i/c*kt; + ]]> + </operator> + <integration_vectors>E_field</integration_vectors> + <dependencies>density_matrix</dependencies> + <![CDATA[ + // Density matrix is Hermitian so we use r_ij=conj(r_ji) + r1301 = conj(r0113); + r1402 = conj(r0214); + r0903 = conj(r0309); + r1503 = conj(r0315); + r1004 = conj(r0410); + r1604 = conj(r0416); + r1105 = conj(r0511); + r0206 = conj(r0602); + r1406 = conj(r0614); + r0307 = conj(r0703); + r0907 = conj(r0709); + r1507 = conj(r0715); + r0408 = conj(r0804); + r1008 = conj(r0810); + r1608 = conj(r0816); + r1509 = conj(r0915); + r1610 = conj(r1016); + + dE1_dz = 0.16666666666666666*i*eta1*(2.449489742783178*r0903 + 4.242640687119286*r1004 + 6.*r1105) - Lt[E1]; + dE2_dz = -0.4082482904638631*i*eta1*(r0907 + r1008) - Lt[E2]; + dE3_dz = -0.3333333333333333*i*eta2*(1.7320508075688772*r1406 + 3.*r1507 + 4.242640687119286*r1608) - Lt[E3]; + dE4_dz = (i*eta2*(2.449489742783178*r1301 + 3*r1402 + 3*r1503 + 2.449489742783178*r1604))/3. - Lt[E4]; + ]]> + </operators> + </integrate> + </sequence> + + + + <!-- The output to generate --> + <output format="binary" filename="realistic_Rb.xsil"> + <group> + <sampling basis="t(1000) " initial_sample="yes"> + <dependencies>E_field_avgd</dependencies> + <moments>I1_out I2_out I3_out I4_out</moments> + <![CDATA[ + I1_out = mod2(E1a); + I2_out = mod2(E2a); + I3_out = mod2(E3a); + I4_out = mod2(E4a); + ]]> + </sampling> + </group> + + <!-- + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>density_matrix_averaged</dependencies> + <moments> + r11_out r22_out r33_out r44_out + r12_re_out r12_im_out r13_re_out r13_im_out r14_re_out r14_im_out + r23_re_out r23_im_out r24_re_out r24_im_out + r34_re_out r34_im_out + </moments> + <![CDATA[ + // populations output + r11_out = r11a.Re(); + r22_out = r22a.Re(); + r33_out = r33a.Re(); + r44_out = r44a.Re(); + // coherences output + r12_re_out = r12a.Re(); + r12_im_out = r12a.Im(); + r13_re_out = r13a.Re(); + r13_im_out = r13a.Im(); + r14_re_out = r14a.Re(); + r14_im_out = r14a.Im(); + r23_re_out = r23a.Re(); + r23_im_out = r23a.Im(); + r24_re_out = r24a.Re(); + r24_im_out = r24a.Im(); + r34_re_out = r34a.Re(); + r34_im_out = r34a.Im(); + ]]> + </sampling> + </group> + --> + + <!-- use the following two groups only for debuging + otherwise they are quite useless and have to much information + in 3D space (z,t,v) --> + <!-- + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>E_field</dependencies> + <moments>I1_out_v I2_out_v I3_out_v I4_out_v</moments> + <![CDATA[ + // light field intensity distribution in velocity subgroups + I1_out_v = mod2(E1); + I2_out_v = mod2(E2); + I3_out_v = mod2(E3); + I4_out_v = mod2(E4); + ]]> + </sampling> + </group> + + <group> + <sampling basis="t(100) v(10)" initial_sample="yes"> + <dependencies>density_matrix</dependencies> + <moments> + r11_out_v r22_out_v r33_out_v r44_out_v + r12_re_out_v r12_im_out_v r13_re_out_v r13_im_out_v r14_re_out_v r14_im_out_v + r23_re_out_v r23_im_out_v r24_re_out_v r24_im_out_v + r34_re_out_v r34_im_out_v + </moments> + <![CDATA[ + // density matrix distribution in velocity subgroups + // populations output + r11_out_v = r11.Re(); + r22_out_v = r22.Re(); + r33_out_v = r33.Re(); + r44_out_v = r44.Re(); + // coherences output + r12_re_out_v = r12.Re(); + r12_im_out_v = r12.Im(); + r13_re_out_v = r13.Re(); + r13_im_out_v = r13.Im(); + r14_re_out_v = r14.Re(); + r14_im_out_v = r14.Im(); + r23_re_out_v = r23.Re(); + r23_im_out_v = r23.Im(); + r24_re_out_v = r24.Re(); + r24_im_out_v = r24.Im(); + r34_re_out_v = r34.Re(); + r34_im_out_v = r34.Im(); + ]]> + </sampling> + </group> + --> + + </output> + + +<info> +Script compiled with XMDS2 version VERSION_PLACEHOLDER (SUBVERSION_REVISION_PLACEHOLDER) +See http://www.xmds.org for more information. + +Variables that can be specified on the command line: + Command line argument E1o = 0.000000e+00 + Command line argument E2o = 1.000000e+02 + Command line argument E3o = 0.000000e+00 + Command line argument E4o = 0.000000e+00 + Command line argument delta1 = 0.000000e+00 + Command line argument delta2 = 0.000000e+00 + Command line argument delta3 = 0.000000e+00 + Command line argument Pwidth = 4.000000e-07 + Command line argument Lcell = 1.000000e-01 + Command line argument Ndens = 1.000000e+15 + Command line argument Temperature = 1.000000e-03 +</info> + +<XSIL Name="moment_group_1"> + <Param Name="n_independent">2</Param> + <Array Name="variables" Type="Text"> + <Dim>6</Dim> + <Stream><Metalink Format="Text" Delimiter=" \n"/> +z t I1_out I2_out I3_out I4_out + </Stream> + </Array> + <Array Name="data" Type="double"> + <Dim>101</Dim> + <Dim>1000</Dim> + <Dim>6</Dim> + <Stream><Metalink Format="Binary" UnsignedLong="uint32" precision="double" Type="Remote" Encoding="LittleEndian"/> +realistic_Rb_expected_mg0.dat + </Stream> + </Array> +</XSIL> +</simulation> diff --git a/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat Binary files differnew file mode 100644 index 0000000..ff30ab5 --- /dev/null +++ b/xmds2/realistic_Rb/tests/testsuite/realistic_Rb_expected_mg0.dat |