diff options
Diffstat (limited to 'Nlevels_with_MOR.xmds')
| -rw-r--r-- | Nlevels_with_MOR.xmds | 481 |
1 files changed, 0 insertions, 481 deletions
diff --git a/Nlevels_with_MOR.xmds b/Nlevels_with_MOR.xmds deleted file mode 100644 index 99d88a0..0000000 --- a/Nlevels_with_MOR.xmds +++ /dev/null @@ -1,481 +0,0 @@ -<?xml version="1.0"?> -<simulation xmds-version="2"> - - <name>Nlevels_with_MOR</name> - - <author>Eugeniy Mikhailov, M. Guidry</author> - <description> - License GPL. - - Solving split 3-level atom - with field propagation along spatial axis Z - with Doppler broadening. - - - * - * -------- |a> - * / / \ \ - * EdL / / \ \ - * / / \ \ - * |c> ----------/--- \ \ EpR - * / EpL \ \ - * |C> -------------- \ \ - * EdR \ \ - * \ \ - * ------\------ |b> - * \ - * ------------- |B> - * - * - - - 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[ - const double pi = M_PI; - 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 - const double Kvec = 2*M_PI/lambda; // k-vector - 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)] - - // --------- 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 Ga=3.0 *(2*M_PI*1e6); - //const double G4=3.0 *(2*M_PI*1e6); - - complex g = Gamma_super; - complex gbc = (2*M_PI)*1e3; // 1/s untits - const complex Split = 0; - const complex noise = 0; - - complex Gab, GAB, Gca, GCA, Gcb, GCB; - - // total decay of i-th level branching ratios. Rij branching of i-th level to j-th - //const double Rab=0.5, Rac=0.5; - - - complex EdLac, EdRac, EpLac, EpRac; - - // 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="EdLo" type="real" default_value="2*1.5*(2*M_PI*1e6)" /> - <argument name="EdRo" type="real" default_value="2*1.5*(2*M_PI*1e6)" /> - - <argument name="EpLo" type="real" default_value="0.05*(2*M_PI*1e6)" /> - <argument name="EpRo" type="real" default_value="0.05*(2*M_PI*1e6)" /> - - - <!-- Fields detuning in [1/s] --> - <argument name="delta_dL" type="real" default_value="0.0" /> - <argument name="delta_pL" type="real" default_value="0.0" /> - <argument name="delta_dR" type="real" default_value="0.0" /> - <argument name="delta_pR" 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; - ]]> - </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>EdL EdR EpL EpR</components> - <initialisation> - <![CDATA[ - // Initial (at starting 'z' position) electromagnetic field does not depend on detuning - // as well as time - EdL=EdLo*exp(-pow( ((t-0.0)/Pwidth),2) ); - EdR=EdRo*exp(-pow( ((t-0.0)/Pwidth),2) ); - EpL=EpLo*exp(-pow( ((t-0.0)/Pwidth),2) ); - EpR=EpRo*exp(-pow( ((t-0.0)/Pwidth),2) ); - - //EpL = EpLo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5)); - //EpR = EpRo*(1+0.01*2*(((double)rand() / (double)RAND_MAX)-0.5)); - - //E2=E2o*exp(-pow( ((t-0.0)/Pwidth),2) ); - ]]> - </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>EdLa EdRa EpLa EpRa</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; - - EdLa=EdL*prob_v_normalized; - EdRa=EdR*prob_v_normalized; - EpLa=EpL*prob_v_normalized; - EpRa=EpR*prob_v_normalized; - ]]> - </evaluation> - </computed_vector> - - <!-- Averaged across Maxwell distribution density matrix components --> - <computed_vector name="density_matrix_averaged" dimensions="t" type="complex"> - <components>rbb_av rBB_av rcc_av rCC_av raa_av rcb_av rab_av rca_av rCB_av rAB_av rCA_av</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; - - rbb_av=rbb*prob_v_normalized; - rBB_av=rBB*prob_v_normalized; - rcc_av=rcc*prob_v_normalized; - rCC_av=rCC*prob_v_normalized; - raa_av=raa*prob_v_normalized; - - rcb_av=rcb*prob_v_normalized; - rab_av=rab*prob_v_normalized; - rca_av=rca*prob_v_normalized; - - rCB_av = rCB*prob_v_normalized; - rAB_av = rAB*prob_v_normalized; - rCA_av = rCA*prob_v_normalized; - - ]]> - </evaluation> - </computed_vector> - - - <vector name="density_matrix" type="complex" initial_space="t"> - <components>rbb rBB rcc rCC raa rcb rab rca rCB rAB rCA</components> - <initialisation> - <!--This sets boundary condition at all times and left border of z (i.e. z=0)--> - <dependencies>E_field_avgd</dependencies> - <![CDATA[ - EdLac = conj(EdLa); - EdRac = conj(EdRa); - EpLac = conj(EpLa); - EpRac = conj(EpRa); - - rbb = 0.25; rcc = 0.25; rBB = 0.25; rCC = 0.25; - raa = 0; - rcb = 0; rab = 0; rca = 0; - rCB = 0; rAB = 0; rCA = 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 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))--> - </boundary_condition> - <![CDATA[ - - Gab=g+i*(delta_dL+delta_pR+0*noise); - GAB=g+i*(delta_pL+delta_dR+0*noise); - - Gca=g-i*(delta_dL); - GCA=g-i*(delta_pL); - - Gcb=gbc+i*(Split + delta_pR+0*noise); - GCB=gbc+i*(-Split + delta_pL+0*noise); - - complex rba=conj(rab); - complex rac=conj(rca); - complex rbc=conj(rcb); - complex rBA=conj(rAB); - complex rAC=conj(rCA); - complex rBC=conj(rCB); - - draa_dt = -i*EpRac*rab+i*EpRa*rba-i*EdLac*rac+i*EdLa*rca-2*g*raa - -i*EdRac*rAB+i*EdRa*rBA-i*EpLac*rAC+i*EpLa*rCA; - - drbb_dt = i*EpRac*rab-i*EpRa*rba+g*raa-gbc*rbb+gbc*rcc; - drBB_dt = i*EdRac*rAB-i*EdRa*rBA+g*raa-gbc*rBB+gbc*rCC; - - drcc_dt = i*EdLac*rac-i*EdLa*rca+g*raa-gbc*rcc+gbc*rbb; - drCC_dt = i*EpLac*rAC-i*EpLa*rCA+g*raa-gbc*rCC+gbc*rBB; - - drab_dt = -Gab*rab+i*EpRa*(rbb-raa)+i*EdLa*rcb; - drca_dt = -Gca*rca+i*EdLac*(raa-rcc)-i*EpRac*rcb; - drcb_dt = -Gcb*rcb-i*EpRa*rca+i*EdLac*rab; - - drAB_dt = -Gab*rAB+i*EdRa*(rBB-raa)+i*EpLa*rCB; - drCA_dt = -Gca*rCA+i*EpLac*(raa-rCC)-i*EdRac*rCB; - drCB_dt = -GCB*rCB-i*EdRa*rCA+i*EpLac*rAB; - - - ]]> - </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[ - dEdL_dz = i*eta*(rca) +0*Lt[EdL] ; - dEdR_dz = i*eta*(rAB) +0*Lt[EdR] ; - dEpL_dz = i*eta*(rCA) +0*Lt[EpL] ; - dEpR_dz = i*eta*(rab) +0*Lt[EpR] ; - - ]]> - </operators> - </integrate> - </sequence> - - - - <!-- The output to generate --> - <output format="binary" filename="Nlevels_with_MOR.xsil"> - <group> - <sampling basis="t(1000) " initial_sample="yes"> - <dependencies>E_field_avgd</dependencies> - <moments>IdL_out IpL_out IdR_out IpR_out</moments> - <![CDATA[ - IdL_out = mod2(EdLa); - IpL_out = mod2(EpLa); - IdR_out = mod2(EdRa); - IpR_out = mod2(EpRa); - ]]> - </sampling> - </group> - - <group> - <sampling basis="t(100) v(2)" initial_sample="yes"> - <dependencies>density_matrix_averaged</dependencies> - <moments> - rbb_out rBB_out - rcc_out rCC_out - raa_out - rcb_re_out rcb_im_out - rCB_re_out rCB_im_out - rab_re_out rab_im_out - rAB_re_out rAB_im_out - rca_re_out rca_im_out - rCA_re_out rCA_im_out - </moments> - <![CDATA[ - // populations output - rbb_out = rbb_av.Re(); - rBB_out = rBB_av.Re(); - - rcc_out = rcc_av.Re(); - rCC_out = rCC_av.Re(); - - raa_out = raa_av.Re(); - - // coherences output - rcb_re_out = rcb_av.Re(); - rcb_im_out = rcb_av.Im(); - - rCB_re_out = rCB_av.Re(); - rCB_im_out = rCB_av.Im(); - - rab_re_out = rab_av.Re(); - rab_im_out = rab_av.Im(); - - rAB_re_out = rAB_av.Re(); - rAB_im_out = rAB_av.Im(); - - rca_re_out = rca_av.Re(); - rca_im_out = rca_av.Im(); - - rCA_re_out = rCA_av.Re(); - rCA_im_out = rCA_av.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> - rbb_out_v rcc_out_v raa_out_v r44_out_v - rbc_re_out_v rbc_im_out_v rba_re_out_v rba_im_out_v r14_re_out_v r14_im_out_v - rca_re_out_v rca_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 - rbb_out_v = rbb.Re(); - rcc_out_v = rcc.Re(); - raa_out_v = raa.Re(); - r44_out_v = r44.Re(); - // coherences output - rbc_re_out_v = rbc.Re(); - rbc_im_out_v = rbc.Im(); - rba_re_out_v = rba.Re(); - rba_im_out_v = rba.Im(); - r14_re_out_v = r14.Re(); - r14_im_out_v = r14.Im(); - rca_re_out_v = rca.Re(); - rca_im_out_v = rca.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: ---> |