diff --git a/CMakeLists.txt b/CMakeLists.txt index e3cf45c..b1e2d73 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,6 +1,6 @@ cmake_minimum_required(VERSION 3.14 FATAL_ERROR) -project(JoSIM VERSION 2.5.8) +project(JoSIM VERSION 2.6) set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) diff --git a/README.md b/README.md index 43257df..3453844 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # JoSIM ### Superconductor Circuit Simulator -##### Testing: v2.5.8 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoSIM-CI-Devel?branchName=testing)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=1&branchName=testing) +##### Testing: v2.6 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoSIM-CI-Devel?branchName=testing)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=1&branchName=testing) -##### Stable: v2.5.8 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoeyDelp.JoSIM?branchName=master)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=3&branchName=master) +##### Stable: v2.6 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoeyDelp.JoSIM?branchName=master)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=3&branchName=master) --- @@ -35,6 +35,11 @@ Referencing: --- ## Changelog +### v2.6 +- Full Johnson-Nyquist temperature noise for Resistors and Resistive branch of RCSJ model +- Phi-junction model with verified Pi junction results +- Various bug fixes and code clean up + ### v2.5.8 - Updated KLU to a 64-bit version capable of handling non-zero count larger than the maximum value of a 32-bit integer. - Included SuperLU 64-bit that can be activated using the *-x 1* switch (default 0 for KLU). diff --git a/docs/img/josim_help.png b/docs/img/josim_help.png index 7daf103..e24e4df 100644 Binary files a/docs/img/josim_help.png and b/docs/img/josim_help.png differ diff --git a/docs/img/josim_jtl_ex.png b/docs/img/josim_jtl_ex.png index 4f47bf6..2875c1a 100644 Binary files a/docs/img/josim_jtl_ex.png and b/docs/img/josim_jtl_ex.png differ diff --git a/docs/index.md b/docs/index.md index d3fe10f..a5f2629 100644 --- a/docs/index.md +++ b/docs/index.md @@ -4,9 +4,9 @@ Developers Manual for v2.5.1 ## Project Status -### Testing: v2.5.1 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoSIM-CI-Devel?branchName=testing)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=1&branchName=testing) +### Testing: v2.6 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoSIM-CI-Devel?branchName=testing)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=1&branchName=testing) -### Stable: v2.5.1 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoeyDelp.JoSIM?branchName=master)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=3&branchName=master) +### Stable: v2.6 - Status: [![Build Status](https://joeydelp.visualstudio.com/JoSIM/_apis/build/status/JoeyDelp.JoSIM?branchName=master)](https://joeydelp.visualstudio.com/JoSIM/_build/latest?definitionId=3&branchName=master) ## Introduction diff --git a/docs/syntax.md b/docs/syntax.md index f53356f..ed4cbb4 100644 --- a/docs/syntax.md +++ b/docs/syntax.md @@ -67,7 +67,7 @@ The value of a capacitor is in Farad. ### Josephson Junction (JJ) -**B**Label $N^{+}$ $N^{-}$ ** **MODEL** [area=<**AREA**\>] [ic=<**IC**>] +**B**Label $N^{+}$ $N^{-}$ ** **MODEL** [area=<**AREA**\>] [ic=<**IC**>] [temp=<**TEMP**>] [neb=<**FREQ**>] A Josephson junction is a two terminal device but could also be defined with a third non-connected node to allow compatibility with WRspice. This node is not used in any way in JoSIM. @@ -75,6 +75,8 @@ The Josephson junction requires specification of a model name which can be defin When **AREA** or **IC** is not specified then an area=1 is used as default. +The temp and neb commands have the same descriptions as for the resistor. + #### Model This model control has the following syntax @@ -96,12 +98,14 @@ The only junction model currently supported by JoSIM is the RCSJ model and thus | DELV | 0, $\infty$ | 0.1E-3 | Transitional voltage from subgap to normal | | D | 0.0, 1 | 0.0 | Point of contact transparency affecting current phase relationship | | ICFCT or ICFACT | 0, 1 | $\frac{\pi}{4}$ | Ratio of critical current to quasiparticle step height | -| PHI | 0, $\pi$ | 0 | Starting phase of junction | +| PHI | 0, $2\pi$ | 0 | Allows phi-junction capability such as the $\pi$-junction. | The *.model* line is unique to the subcircuit it falls under and can thus allow different models with the same name under separate subcircuits. If the model is not found under the subcircuit it will be searched for globally and if not found default values (default model) will be used instead. The **AREA** and **IC** parameters act as modifiers to the model parameters. **AREA** is a critical current multiplier, where if **IC** is specified it replaces the **AREA** value by $AREA=\frac{IC_{jj}}{IC_{model}}$. +By setting the **PHI** parameter of the model, the phase value is persistantly subtracted from the phase ($\phi$) in the $\sin(\phi)$ part of the JJ current. This allows elements such as the $\pi$-junction to be modeled. + ### Transmission Line **T**Label $N^{+}_{1}$ $N^{-}_{1}$ $N^{+}_{2}$ $N^{-}_{2}$ **TD=VALUE** **Z0=VALUE** @@ -247,12 +251,14 @@ The most important of these control commands is the transient simulation command ### Transient Analysis -**.tran** $T_{STEP}$ $T_{STOP}$ [$P_{START}$ [$P_{STEP}$]] +**.tran** $T_{STEP}$ $T_{STOP}$ [$P_{START}$ [$P_{STEP}$]] DST This generates a simulation that runs from 0 until $T_{STOP}$. The amount simulation steps that will be performed is $n=\frac{T_{STOP}}{T_{STEP}}$. $P_{START}$ indicates at what point output will start printing. $P_{STEP}$ sets the size of the print steps. This has to be larger or equal to $T_{STEP}$. +DST disables the start-up time. The start-up time is a period calculated internally by the simulator in which components settle. This is equivalent to the few picoseconds from when a circuit initially receives power (power switch flipped). + ### Subcircuits Subcircuits allow subdivision and reuse of smaller circuits within a larger design. When wrapped in a subcircuit control devices are allowed to have the same label names as specified elsewhere in the netlist as the subcircuit completely isolates them. @@ -273,6 +279,8 @@ A subcircuit can be used in the main netlist or another subcircuit (nesting) usi **X**Label *IO Nodes* *SubcktName* (WRspice (normal SPICE) mode) +Additionally, keywords in the form of **LABEL=VALUE** can be appended to the end of the subcircuit declaration line which when instantiated will replace the value of the **LABEL** component within the subcircuit with the associated **VALUE**. This allows for unique subcircuit instantiations which would prove useful in testing various parameters without altering the original subcircuit or having multiple instances of the same subcircuit definition. This could open the door for potential future margin and optimization software. + ### Noise As mentioned in the technical discussion, noise can be automatically inserted as current sources in parallel to each resistor. This thermal noise temperature and bandwidth can be specified globally using the following commands: diff --git a/docs/tech_disc.md b/docs/tech_disc.md index 09089de..1f58cda 100644 --- a/docs/tech_disc.md +++ b/docs/tech_disc.md @@ -58,6 +58,10 @@ This enables different levels of verbosity by the simulator. Level `1` verbosity Only displays the version string for **josim-cli** and then exits. This string is always displayed by default for any command. This command, however, only displays the version string. +### Solver (-x): + +This option allows switching between LU solvers. Default is KLU solver. Setting this to `1` enables SuperLU solver. + ## Parse Input {#parse-input} The input data received through any of the two means above is parsed in the following steps: diff --git a/include/JoSIM/BasicComponent.hpp b/include/JoSIM/BasicComponent.hpp index 08f0647..37ac01d 100644 --- a/include/JoSIM/BasicComponent.hpp +++ b/include/JoSIM/BasicComponent.hpp @@ -27,7 +27,7 @@ namespace JoSIM { class IndexInfo { public: int_o posIndex_, negIndex_, currentIndex_; - NodeConfig nodeConfig_; + NodeConfig nodeConfig_ = NodeConfig::GND; }; class MatrixInfo { diff --git a/include/JoSIM/Errors.hpp b/include/JoSIM/Errors.hpp index 6394502..cf6bf0c 100644 --- a/include/JoSIM/Errors.hpp +++ b/include/JoSIM/Errors.hpp @@ -35,6 +35,7 @@ namespace JoSIM { MISSING_MAIN, UNKNOWN_SUBCKT, EMPTY_FILE, + IO_MISMATCH, UNKNOWN_CONTROL }; @@ -188,10 +189,10 @@ namespace JoSIM { static void control_errors( ControlErrors errorCode, string_o message = std::nullopt); - [[noreturn]] static void model_errors( + static void model_errors( ModelErrors errorCode, string_o message = std::nullopt); - static void matrix_errors( + [[noreturn]] static void matrix_errors( MatrixErrors errorCode, string_o message = std::nullopt); [[noreturn]] static void misc_errors( @@ -209,7 +210,7 @@ namespace JoSIM { static void output_errors( OutputErrors errorCode, string_o message = std::nullopt); - static void netlist_errors( + [[noreturn]] static void netlist_errors( NetlistErrors errorCode, string_o message = std::nullopt); static void verbosity_errors( diff --git a/include/JoSIM/Function.hpp b/include/JoSIM/Function.hpp index e9c8c84..5527487 100644 --- a/include/JoSIM/Function.hpp +++ b/include/JoSIM/Function.hpp @@ -25,6 +25,7 @@ namespace JoSIM { FunctionType fType_ = FunctionType::PWL; std::vector timeValues_; std::vector ampValues_; + std::vector miscValues_; void parse_pwl(const tokens_t& t, const Input& iObj, const string_o& s); void parse_pulse(const tokens_t& t, const Input& iObj, const string_o& s); void parse_sin(const tokens_t& t, const Input& iObj, const string_o& s); @@ -46,6 +47,12 @@ namespace JoSIM { const string_o& subckt); double value(double x); void ampValues(std::vector values); + std::vector ampValues() { + return ampValues_; + } + void clearMisc() { + miscValues_.clear(); + } }; // class Function diff --git a/include/JoSIM/Input.hpp b/include/JoSIM/Input.hpp index 525b79b..de7d2c3 100644 --- a/include/JoSIM/Input.hpp +++ b/include/JoSIM/Input.hpp @@ -24,7 +24,8 @@ namespace JoSIM { public: Netlist netlist; Transient transSim; - std::optional globalTemp, neB; + std::optional globalTemp; + double neB = 1E12; std::vector fileLines, controls; std::vector relevantX; std::unordered_map parameters; @@ -32,8 +33,8 @@ namespace JoSIM { std::vector output_files; std::optional fileParentPath; - Input() {}; - Input(CliOptions& cli_options) { + Input() : netlist(*this) {}; + Input(CliOptions& cli_options) : netlist(*this) { argAnal = cli_options.analysis_type; argVerb = cli_options.verbose; argMin = cli_options.minimal; diff --git a/include/JoSIM/JJ.hpp b/include/JoSIM/JJ.hpp index c5325ac..02a150f 100644 --- a/include/JoSIM/JJ.hpp +++ b/include/JoSIM/JJ.hpp @@ -10,6 +10,7 @@ #include "JoSIM/AnalysisType.hpp" #include "JoSIM/Input.hpp" #include "JoSIM/Spread.hpp" +#include "JoSIM/Function.hpp" #include #include @@ -54,16 +55,17 @@ namespace JoSIM { public: int variableIndex_ = 0; double area_ = 1.0; - std::optional Ic_; - std::optional model_; + std::optional Ic_, temp_, neb_, spAmp_; + Model model_; double phaseConst_ = 0.0; double lowerB_ = 0.0, upperB_ = 0.0, gLarge_ = 0.0; - double del0_ = 0.0, del_ = 0.0, rncalc_ = 0.0; + double del0_ = 0.0, del_ = 0.0; double pn1_ = 0.0, pn2_ = pn1_, pn3_ = pn2_, pn4_ = pn3_, phi0_ = 0.0; double vn1_ = 0.0, vn2_ = vn1_, vn3_ = vn2_, vn4_ = vn3_, vn5_ = vn4_, vn6_ = vn5_; - double transitionCurrent_ = 0.0; + double it_ = 0.0; JoSIM::AnalysisType at_; + std::optional thermalNoise; JJ( const std::pair& s, const NodeConfig& ncon, @@ -82,8 +84,6 @@ namespace JoSIM { bool update_value(const double& v); - void update_timestep(const double& factor) override; - void step_back() override { pn2_ = pn4_; vn3_ = vn6_; diff --git a/include/JoSIM/LUSolve.hpp b/include/JoSIM/LUSolve.hpp index ec8ec55..0ff9b56 100644 --- a/include/JoSIM/LUSolve.hpp +++ b/include/JoSIM/LUSolve.hpp @@ -12,7 +12,6 @@ namespace JoSIM { class LUSolve { private: char equed[1] = { 'N' }; - yes_no_t equil; trans_t trans; SuperMatrix A, L, U, B, X; GlobalLU_t Glu; @@ -21,11 +20,11 @@ namespace JoSIM { long long* etree; void* work = nullptr; long long info, lwork, nrhs, ldx; - long long i, m, n, nnz; + long long m, n, nnz; double* rhsb,* rhsx,* xact; double* R,* C; double* ferr,* berr; - double u, rpg, rcond; + double rpg, rcond; mem_usage_t mem_usage; superlu_options_t options; SuperLUStat_t stat; diff --git a/include/JoSIM/Matrix.hpp b/include/JoSIM/Matrix.hpp index b40dacc..9e40c1d 100644 --- a/include/JoSIM/Matrix.hpp +++ b/include/JoSIM/Matrix.hpp @@ -16,6 +16,9 @@ namespace JoSIM { class Matrix { + private: + std::vector nodeConfig, nodeConfig2; + bool needsTR_ = true; public: AnalysisType analysisType = AnalysisType::Phase; std::vector sourcegen; @@ -25,19 +28,22 @@ namespace JoSIM { std::unordered_map nm; nodeconnections nc; std::unordered_set lm; - int branchIndex; - std::vector nz, nz_orig; + int branchIndex = 0; + std::vector nz; std::vector ci, rp; std::vector relevantTraces; std::vector relevantIndices; Matrix() {}; void create_matrix(Input& iObj); + void setup(Input& iObj); + void create_components(Input& iObj); + void handle_mutual_inductance(Input& iObj); + void reduce_step(Input& iObj); void create_csr(); void create_nz(); void create_ci(); void create_rp(); - void mod_timestep(const int factor); }; } // namespace JoSIM #endif diff --git a/include/JoSIM/Misc.hpp b/include/JoSIM/Misc.hpp index b37ef7c..288ca12 100644 --- a/include/JoSIM/Misc.hpp +++ b/include/JoSIM/Misc.hpp @@ -78,7 +78,9 @@ namespace JoSIM { template std::string precise_to_string(const T a_value) { std::ostringstream out; - out << std::uppercase << std::fixed << std::scientific << a_value; + out << std::uppercase << + std::setprecision(std::numeric_limits::digits10 + 1) << + a_value; return out.str(); } diff --git a/include/JoSIM/Model.hpp b/include/JoSIM/Model.hpp index d07565b..5a2fb90 100644 --- a/include/JoSIM/Model.hpp +++ b/include/JoSIM/Model.hpp @@ -11,114 +11,120 @@ namespace JoSIM { class Model { private: std::string modelName_; - double voltageGap_; - double criticalCurrent_; - int resistanceType_; - double normalResistance_; - double subgapResistance_; - double capacitance_; - double temperature_; - double criticalTemperature_; + double vg_; + double ic_; + int rtype_; + double rn_; + double r0_; + double c_; + double t_; + double tc_; double deltaV_; - double transparency_; - double criticalToNormalRatio_; - double phaseOffset_; + double d_; + double icFct_; + double phiOff_; + bool tDep_; public: Model() : - voltageGap_(2.8E-3), - criticalCurrent_(1E-3), - resistanceType_(1), - normalResistance_(5), - subgapResistance_(30), - capacitance_(2.5E-12), - temperature_(4.2), - criticalTemperature_(9.1), + vg_(2.8E-3), + ic_(1E-3), + rtype_(1), + rn_(5), + r0_(30), + c_(2.5E-12), + t_(4.2), + tc_(9.1), deltaV_(0.1E-3), - transparency_(0), - criticalToNormalRatio_(Constants::PI / 4), - phaseOffset_(0) {}; + d_(0), + icFct_(Constants::PI / 4), + phiOff_(0), + tDep_(false) {}; - std::string get_modelName() const { + std::string modelName() const { return modelName_; } - double get_voltageGap() const { - return voltageGap_; + void modelName(const std::string& n) { + modelName_ = n; } - double get_criticalCurrent() const { - return criticalCurrent_; + double vg() const { + return vg_; } - int get_resistanceType() const { - return resistanceType_; + void vg(const double& v) { + vg_ = v; } - double get_normalResistance() const { - return normalResistance_; + double ic() const { + return ic_; } - double get_subgapResistance() const { - return subgapResistance_; + void ic(const double& i) { + ic_ = i; } - double get_capacitance() const { - return capacitance_; + int rtype() const { + return rtype_; } - double get_temperature() const { - return temperature_; + void rtype(const int& r) { + rtype_ = r; } - double get_criticalTemperature() const { - return criticalTemperature_; + double rn() const { + return rn_; } - double get_deltaV() const { - return deltaV_; + void rn(const double& r) { + rn_ = r; } - double get_transparency() const { - return transparency_; + double r0() const { + return r0_; } - double get_criticalToNormalRatio() const { - return criticalToNormalRatio_; + void r0(const double& r) { + r0_ = r; } - double get_phaseOffset() const { - return phaseOffset_; + double c() const { + return c_; } - - void set_modelName(const std::string& n) { - modelName_ = n; + void c(const double& c) { + c_ = c; + } + double t() const { + return t_; } - void set_voltageGap(const double& v) { - voltageGap_ = v; + void t(const double& t) { + t_ = t; } - void set_criticalCurrent(const double& i) { - criticalCurrent_ = i; + double tc() const { + return tc_; } - void set_resistanceType(const int& r) { - resistanceType_ = r; + void tc(const double& t) { + tc_ = t; } - void set_normalResistance(const double& r) { - normalResistance_ = r; + double deltaV() const { + return deltaV_; } - void set_subgapResistance(const double& r) { - subgapResistance_ = r; + void deltaV(const double& d) { + deltaV_ = d; } - void set_capacitance(const double& c) { - capacitance_ = c; + double d() const { + return d_; } - void set_temperature(const double& t) { - temperature_ = t; + void d(const double& t) { + d_ = t; } - void set_criticalTemperature(const double& t) { - criticalTemperature_ = t; + double icFct() const { + return icFct_; } - void set_deltaV(const double& d) { - deltaV_ = d; + void icFct(const double& r) { + icFct_ = r; } - void set_transparency(const double& t) { - transparency_ = t; + double phiOff() const { + return phiOff_; } - void set_criticalToNormalRatio(const double& r) { - criticalToNormalRatio_ = r; + void phiOff(const double& o) { + phiOff_ = o; } - void set_phaseOffset(const double& o) { - phaseOffset_ = o; + bool tDep() { + return tDep_; + } + void tDep(bool b) { + tDep_ = b; } - static void parse_model( const std::pair& s, vector_pair_t& models, const param_map& p); diff --git a/include/JoSIM/Netlist.hpp b/include/JoSIM/Netlist.hpp index 68bb96b..64efda7 100644 --- a/include/JoSIM/Netlist.hpp +++ b/include/JoSIM/Netlist.hpp @@ -12,6 +12,8 @@ #include #include +typedef std::unordered_map s_map; + namespace JoSIM { struct pair_hash { template @@ -38,15 +40,20 @@ namespace JoSIM { }; }; + class Input; + class Netlist { + private: + Input& parent; void id_io_subc_label( - const tokens_t& lineTokens, tokens_t& io, + const tokens_t& lineTokens, tokens_t& io, s_map& params, std::string& subcktName, std::string& label, const std::unordered_map& subcircuits); bool rename_io_nodes( std::string& node, const tokens_t& subIO, const tokens_t& parentIO); - void expand_io( - Subcircuit& subc, const tokens_t& io, const std::string& label); + void expand_io(Subcircuit& subc, s_map& params, const tokens_t& io, + const std::string& label); + void insert_parameter(tokens_t& t, s_map& params); public: std::unordered_map< std::pair, tokens_t, pair_hash> models; @@ -58,7 +65,8 @@ namespace JoSIM { std::vector> expNetlist; int jjCount, compCount, subcktCounter, nestedSubcktCount, subcktTotal = 0; bool containsSubckt, argMin = false; - Netlist() : + Netlist(Input& p) : + parent(p), jjCount(0), compCount(0), subcktCounter(0), diff --git a/include/JoSIM/Noise.hpp b/include/JoSIM/Noise.hpp index b58813c..f626cef 100644 --- a/include/JoSIM/Noise.hpp +++ b/include/JoSIM/Noise.hpp @@ -8,19 +8,12 @@ namespace JoSIM { class Noise { - private: - static double determine_spectral_amplitude(const double& R, - const double& T); - static std::pair - create_resistive_current_noise(Input& iObj, const double& R, - const double& T, const double& B, - std::vector>::iterator& i); + public: static void determine_global_temperature(Input& iObj); static void determine_noise_effective_bandwidth(Input& iObj); + static double determine_spectral_amplitude(const double& R, + const double& T); - - public: - static void add_noise_sources(Input& iObj); }; } diff --git a/include/JoSIM/Resistor.hpp b/include/JoSIM/Resistor.hpp index cf52c88..b295a8a 100644 --- a/include/JoSIM/Resistor.hpp +++ b/include/JoSIM/Resistor.hpp @@ -9,6 +9,7 @@ #include "JoSIM/AnalysisType.hpp" #include "JoSIM/Input.hpp" #include "JoSIM/Spread.hpp" +#include "JoSIM/Function.hpp" #include #include @@ -36,8 +37,10 @@ namespace JoSIM { class Resistor : public BasicComponent { private: JoSIM::AnalysisType at_; + std::optional spAmp_, temp_, neb_; public: double pn1_ = 0.0, pn2_ = 0.0, pn3_ = 0.0, pn4_ = 0.0; + std::optional thermalNoise; Resistor( const std::pair& s, const NodeConfig& ncon, diff --git a/include/JoSIM/Simulation.hpp b/include/JoSIM/Simulation.hpp index 645bd96..e6ccb75 100644 --- a/include/JoSIM/Simulation.hpp +++ b/include/JoSIM/Simulation.hpp @@ -30,12 +30,13 @@ namespace JoSIM { class Simulation { private: bool SLU = false; - std::vector x_, b_, xn2_, xn3_; + std::vector x_, b_; int simSize_; JoSIM::AnalysisType atyp_; bool minOut_; bool needsLU_; - bool needsTR_; + bool needsTR_ = true; + bool startup_; double stepSize_, prstep_, prstart_; int simOK_; klu_l_symbolic* Symbolic_; @@ -43,13 +44,13 @@ namespace JoSIM { klu_l_numeric* Numeric_; LUSolve lu; + void setup(Input& iObj, Matrix& mObj); void trans_sim(Matrix& mObj); void setup_b(Matrix& mObj, int i, double step, double factor = 1); - void reduce_step(Matrix& mObj, double factor, - int& stepCount, double currentStep); + void reduce_step(Input& iObj, Matrix& mObj); void handle_cs(Matrix& mObj, double& step, const int& i); - void handle_resistors(Matrix& mObj); + void handle_resistors(Matrix& mObj, double& step); void handle_inductors(Matrix& mObj, double factor = 1); void handle_capacitors(Matrix& mObj); void handle_jj(Matrix& mObj, int& i, double& step, double factor = 1); diff --git a/include/JoSIM/Transient.hpp b/include/JoSIM/Transient.hpp index dc43dfb..bd507cf 100644 --- a/include/JoSIM/Transient.hpp +++ b/include/JoSIM/Transient.hpp @@ -15,12 +15,14 @@ namespace JoSIM { double tstop_; double prstart_; double prstep_; + bool startup_; public: Transient() : tstep_(0.25E-12), tstop_(0.0), prstart_(0.0), - prstep_(1E-12) {}; + prstep_(1E-12), + startup_(true) {}; double tstep() const { return tstep_; @@ -37,6 +39,9 @@ namespace JoSIM { int simsize() const { return static_cast(tstop_ / tstep_); } + bool startup() const { + return startup_; + } void tstep(double value) { tstep_ = value; @@ -50,9 +55,12 @@ namespace JoSIM { void prstep(double value) { prstep_ = value; } + void startup(bool value) { + startup_ = value; + } static void identify_simulation( - const std::vector& controls, Transient& tObj); + std::vector& controls, Transient& tObj); }; } diff --git a/include/JoSIM/TransmissionLine.hpp b/include/JoSIM/TransmissionLine.hpp index 42b1f12..e5c67e0 100644 --- a/include/JoSIM/TransmissionLine.hpp +++ b/include/JoSIM/TransmissionLine.hpp @@ -60,7 +60,7 @@ namespace JoSIM { double n1_1_ = 0.0, n1_2_ = 0.0, n2_1_ = 0.0, n2_2_ = 0.0; double nk_1_ = 0.0, nk_2_ = 0.0, nk1_1_ = 0.0, nk1_2_ = 0.0, nk2_1_ = 0.0, nk2_2_ = 0.0; - int timestepDelay_ = 0.0; + int timestepDelay_ = 0; TransmissionLine( const std::pair& s, const NodeConfig& ncon, diff --git a/scripts/josim-plot b/scripts/josim-plot.py similarity index 67% rename from scripts/josim-plot rename to scripts/josim-plot.py index 4b3b054..b22f5b9 100644 --- a/scripts/josim-plot +++ b/scripts/josim-plot.py @@ -18,7 +18,7 @@ def main(): parser.add_argument("input", help="the CSV input file") parser.add_argument("-o", "--output", help="the output file name with supported extensions: png, jpeg, webp, svg, eps, pdf") parser.add_argument("-x", "--html", help="save the output as an html file for later viewing") - parser.add_argument("-t", "--type", help="type of plot: grid, stacked, combined, square. Default: grid", default="grid") + parser.add_argument("-t", "--type", help="type of plot: grid, stacked, combined, square, sep_comb. Default: grid", default="grid") parser.add_argument("-s", "--subset", nargs='+', help="subset of traces to plot. specify list of column names (as shown in csv file header), ie. \"V(1)\" \"V(2)\". Default: None") parser.add_argument("-c", "--color", help="set the output plot color scheme to one of the following: light, dark, presentation. Default: dark", default='dark') parser.add_argument("-w", "--title", help="set plot title to the provided string") @@ -27,17 +27,20 @@ def main(): # Read arguments from the command line args = parser.parse_args() + # List of possible output formats outformats = [".png", ".jpeg", ".webp", ".svg", ".eps", ".pdf"] - if (os.path.splitext(args.input)[1] == ".csv"): + # How to handle column seperation. CSV uses comma, DAT uses space. + if (os.path.splitext(args.input)[1].lower() == ".csv"): df = pd.read_csv(args.input, sep=',') - elif (os.path.splitext(args.input)[1] == ".dat"): + elif (os.path.splitext(args.input)[1].lower() == ".dat"): df = pd.read_csv(args.input, delim_whitespace=True) else: print("Invalid input file specified: " + args.input) print("Please provide either .csv (comma seperated) or .dat (space seperated) file") sys.exit() + # Determine the plot layout. if(args.type == "grid"): fig = grid_layout(df, args.subset) elif(args.type == "stacked"): @@ -46,11 +49,14 @@ def main(): fig = combined_layout(df, args.subset) elif(args.type == "square"): fig = square_layout(df, args.subset) + elif(args.type == "sep_comb"): + fig = seperate_combined_layout(df, args.subset) else: print("Invalid plot type specified: " + args.type) print("Please provide either gird, stacked or combined as type") sys.exit() - + + # Determine the theme to plot with. if(args.color == 'light'): template = 'plotly_white' elif(args.color == 'dark'): @@ -62,21 +68,37 @@ def main(): print("Please provide either light, dark or presentation as color theme") sys.exit() + # Set the title of the plot if(args.title == None): title=os.path.splitext(os.path.basename(args.input))[0] else: title=args.title + # Update the layout based on the settings fig.update_layout( title=title, title_font_size=30, template=template ) + + # Set the mode bar buttons config = dict({ 'scrollZoom': True, - 'displaylogo': False + 'displaylogo': False, + 'toImageButtonOptions': { + 'format': 'svg', + 'filename': title, + 'height': None, + 'width': None + }, + 'modeBarButtonsToAdd': [ + 'toggleSpikelines', + 'hovercompare', + 'v1hovermode' + ] }) + # Determine wether to show the plot or just dump to file if(args.output == None and args.html == None): fig.show(config=config) elif(args.html != None and args.output == None): @@ -90,11 +112,11 @@ def main(): # Function that sets the Y-axis title relevant to the data def y_axis_title(figLabel): if figLabel[0] == 'V': - return "Voltage (volts)" + return "Voltage (V)" elif figLabel[0] == 'I': - return "Current (ampere)" + return "Current (A)" elif figLabel[0] == 'P': - return "Phase (radians)" + return "Phase (rad)" else: return "Unknown" @@ -191,6 +213,84 @@ def stacked_layout(df, subset): fig['layout']['yaxis' + str(i+1)]['title']=y_axis_title(plots[i]) return fig +# Seperate and Combine like plots +def seperate_combined_layout(df, subset): + plots = df.columns[1:].tolist() if subset == None else subset + V = [] + P = [] + I = [] + U = [] + for i in range(0, len(plots)): + if plots[i][0] == 'V': + V.append(i) + elif plots[i][0] == 'I': + I.append(i) + elif plots[i][0] == 'P': + P.append(i) + else: + U.append(i) + fig_count = 0 + if len(V) != 0: + fig_count = fig_count + 1; + if len(P) != 0: + fig_count = fig_count + 1; + if len(I) != 0: + fig_count = fig_count + 1; + if len(U) != 0: + fig_count = fig_count + 1; + fig = make_subplots( + rows=fig_count, + cols=1, + vertical_spacing=0.2/math.ceil(fig_count/2), + x_title= 'Time (seconds)') + # Add the traces + fig_count = 0 + if len(U) != 0: + fig_count = fig_count + 1; + for i in U: + fig.add_trace(go.Scatter( + x=df.iloc[:,0], y=df.loc[:,plots[i]], + mode='lines', + name=plots[i]), + row=fig_count, + col=1 + ) + fig['layout']['yaxis' + str(fig_count)]['title']=y_axis_title(plots[i]) + if len(V) != 0: + fig_count = fig_count + 1; + for i in V: + fig.add_trace(go.Scatter( + x=df.iloc[:,0], y=df.loc[:,plots[i]], + mode='lines', + name=plots[i]), + row=fig_count, + col=1 + ) + fig['layout']['yaxis' + str(fig_count)]['title']=y_axis_title(plots[i]) + if len(P) != 0: + fig_count = fig_count + 1; + for i in P: + fig.add_trace(go.Scatter( + x=df.iloc[:,0], y=df.loc[:,plots[i]], + mode='lines', + name=plots[i]), + row=fig_count, + col=1 + ) + fig['layout']['yaxis' + str(fig_count)]['title']=y_axis_title(plots[i]) + if len(I) != 0: + fig_count = fig_count + 1; + for i in I: + fig.add_trace(go.Scatter( + x=df.iloc[:,0], y=df.loc[:,plots[i]], + mode='lines', + name=plots[i]), + row=fig_count, + col=1 + ) + fig['layout']['yaxis' + str(fig_count)]['title']=y_axis_title(plots[i]) + return fig + # Combine all the plots def combined_layout(df, subset): plots = df.columns[1:].tolist() if subset == None else subset diff --git a/src/Errors.cpp b/src/Errors.cpp index 56e1b06..ff99530 100644 --- a/src/Errors.cpp +++ b/src/Errors.cpp @@ -157,6 +157,14 @@ void Errors::input_errors(InputErrors errorCode, string_o message) { formattedMessage += "Please check the input file and ensure that the file is not empty."; throw std::runtime_error(formattedMessage); + case InputErrors::IO_MISMATCH: + formattedMessage += + "The IO of line \"" + message.value_or("") + "\" does not " + "match the subcircuit IO.\n"; + formattedMessage += + "Please check the line and ensure correct IO and " + "that parameters do not contain spaces."; + throw std::runtime_error(formattedMessage); case InputErrors::UNKNOWN_CONTROL: formattedMessage += "The control \"" + message.value_or("") + "\" is not known.\n"; @@ -472,7 +480,7 @@ void Errors::control_errors(ControlErrors errorCode, string_o message) { } } -[[noreturn]] void Errors::model_errors(ModelErrors errorCode, +void Errors::model_errors(ModelErrors errorCode, string_o message) { std::string formattedMessage = "Model\n"; switch (errorCode) { @@ -480,8 +488,10 @@ void Errors::control_errors(ControlErrors errorCode, string_o message) { formattedMessage += "Unknown model parameter specified.\n"; formattedMessage += "Model line: " + message.value_or("") + "\n"; formattedMessage += - "Please refer to the model definition in the documentation"; - throw std::runtime_error(formattedMessage); + "Continuing with default model parameters.\n" + "Please consult the syntax guide for more information."; + warning_message(formattedMessage); + break; case ModelErrors::UNKNOWN_MODEL_TYPE: formattedMessage += "Unknown model type specified.\n"; formattedMessage += "Model line: " + message.value_or(""); @@ -497,7 +507,7 @@ void Errors::control_errors(ControlErrors errorCode, string_o message) { } } -void Errors::matrix_errors(MatrixErrors errorCode, string_o message) { +[[noreturn]] void Errors::matrix_errors(MatrixErrors errorCode, string_o message) { std::string formattedMessage = "Matrix\n"; switch (errorCode) { case MatrixErrors::NON_SQUARE: @@ -780,7 +790,7 @@ void Errors::parsing_errors(ParsingErrors errorCode, string_o message) { } } -void Errors::netlist_errors(NetlistErrors errorCode, string_o message) { +[[noreturn]] void Errors::netlist_errors(NetlistErrors errorCode, string_o message) { std::string formattedMessage = "Netlist\n"; switch (errorCode) { case NetlistErrors::NO_SUCH_NODE: @@ -793,6 +803,11 @@ void Errors::netlist_errors(NetlistErrors errorCode, string_o message) { formattedMessage += message.value_or("") + "\n"; formattedMessage += "Please check for any disconnections in the netlist"; throw std::runtime_error(formattedMessage); + default: + formattedMessage += + "Unknown netlist error: " + message.value_or("") + "\n"; + formattedMessage += "Please contact the developer."; + throw std::runtime_error(formattedMessage); } } diff --git a/src/Function.cpp b/src/Function.cpp index 4469307..0daa3e7 100644 --- a/src/Function.cpp +++ b/src/Function.cpp @@ -65,6 +65,8 @@ void Function::parse_function( Misc::tokenize(str).back(), iObj.parameters, subckt))) { fType_ = FunctionType::DC; parse_dc(t, iObj, subckt); + } else { + // Assume DC } } @@ -234,6 +236,7 @@ void Function::parse_noise( const double& tstep = iObj.transSim.tstep(); // Set amplitudes ampValues_.emplace_back(parse_param(t.at(0), iObj.parameters, s)); + ampValues_.emplace_back(0.0); // Default values for the arguments timeValues_.emplace_back(0.0); timeValues_.emplace_back(tstep); @@ -242,6 +245,11 @@ void Function::parse_noise( timeValues_.at(0) = parse_param(t.at(1), iObj.parameters, s); if (t.size() >= 3) timeValues_.at(1) = parse_param(t.at(2), iObj.parameters, s); + miscValues_.emplace_back(0.0); + miscValues_.emplace_back( + ampValues_.at(0) * Misc::grand() / sqrt(2.0 * timeValues_.back())); + miscValues_.emplace_back(timeValues_.at(0)); + miscValues_.emplace_back(timeValues_.at(0) + timeValues_.at(1)); if (iObj.argVerb) { if (ampValues_.at(0) == 0.0) { Errors::function_errors(FunctionErrors::NOISE_VA_ZERO, t.at(1)); @@ -383,10 +391,23 @@ double Function::return_cus(double& x) { } double Function::return_noise(double& x) { - if (x < timeValues_.front()) - return 0.0; - else - return ampValues_.at(0) * Misc::grand() / sqrt(2.0 * timeValues_.back()); + if (x < timeValues_.front()) { + return ampValues_.at(1); + } else { + if (x >= miscValues_.at(3)) { + miscValues_.at(0) = miscValues_.at(1); + miscValues_.at(1) = ampValues_.at(0) * + Misc::grand() / sqrt(2.0 * timeValues_.back()); + miscValues_.at(2) = miscValues_.at(3); + miscValues_.at(3) = miscValues_.at(2) + timeValues_.back(); + } + double& y2 = miscValues_.at(1); + double& y1 = miscValues_.at(0); + double& x2 = miscValues_.at(3); + double& x1 = miscValues_.at(2); + // Calculate function value and return it + return y1 + ((y2 - y1) / (x2 - x1)) * (x - x1); + } } double Function::return_pws(double& x) { diff --git a/src/IV.cpp b/src/IV.cpp index 12af488..21fe9a6 100644 --- a/src/IV.cpp +++ b/src/IV.cpp @@ -41,6 +41,7 @@ void IV::setup_iv(const tokens_t& i, const Input& iObj) { double maxC = parse_param(i.at(2), iObj.parameters); // Create an input object for this simulation Input ivInp = iObj; + ivInp.controls.clear(); tokens_t jj = { "B01", "1", "0", model, "AREA=1" }; tokens_t ib = { "IB01", "0", "1", "PWL(0", "0", "10P", "0", "50P", "2.5U)" }; ivInp.netlist.expNetlist = { @@ -74,24 +75,28 @@ std::vector> IV::generate_iv(double maxC, currentCurr += currentIncrement; ivMat.sourcegen.back().ampValues({ 0, 0, currentCurr }); iv_data.emplace_back(do_simulate(ivInp, ivMat)); + ivInp.transSim.tstep(0.05E-12); } while (currentCurr >= 0) { // Create a matrix object using the input object currentCurr -= currentIncrement; ivMat.sourcegen.back().ampValues({ 0, maxC, currentCurr }); iv_data.emplace_back(do_simulate(ivInp, ivMat)); + ivInp.transSim.tstep(0.05E-12); } while (currentCurr >= -maxC) { // Create a matrix object using the input object currentCurr -= currentIncrement; ivMat.sourcegen.back().ampValues({ 0, 0, currentCurr }); iv_data.emplace_back(do_simulate(ivInp, ivMat)); + ivInp.transSim.tstep(0.05E-12); } while (currentCurr <= 0) { // Create a matrix object using the input object currentCurr += currentIncrement; ivMat.sourcegen.back().ampValues({ 0, -maxC, currentCurr }); iv_data.emplace_back(do_simulate(ivInp, ivMat)); + ivInp.transSim.tstep(0.05E-12); } return iv_data; } diff --git a/src/Input.cpp b/src/Input.cpp index fbb83c6..c47bbcf 100644 --- a/src/Input.cpp +++ b/src/Input.cpp @@ -10,6 +10,7 @@ #include #include #include +#include using namespace JoSIM; @@ -123,6 +124,7 @@ std::vector Input::read_input( } void Input::parse_input(string_o fileName) { + srand(time(NULL)); // Create a seperate thread that will be used for printing creation progress std::thread printingThread; // Variable to store the tokenized input diff --git a/src/JJ.cpp b/src/JJ.cpp index af8882b..ca8ee6a 100644 --- a/src/JJ.cpp +++ b/src/JJ.cpp @@ -5,6 +5,7 @@ #include "JoSIM/Misc.hpp" #include "JoSIM/Errors.hpp" #include "JoSIM/Constants.hpp" +#include "JoSIM/Noise.hpp" #include @@ -42,12 +43,31 @@ JJ::JJ( const nodemap& nm, std::unordered_set& lm, nodeconnections& nc, Input& iObj, Spread& spread, int& bi) { double spr = 1.0; - for (auto i : s.first) { - if (i.find("SPREAD=") != std::string::npos) { - spr = - parse_param(i.substr(i.find("SPREAD=") + 7), iObj.parameters, s.second); + tokens_t t; + for (int i = 3; i < s.first.size(); ++i) { + auto& ti = s.first.at(i); + if (ti.rfind("SPREAD=", 0) == 0) { + spr = parse_param(ti.substr(7), iObj.parameters, s.second); + } else if (ti.rfind("TEMP=", 0) == 0) { + temp_ = parse_param(ti.substr(5), iObj.parameters, s.second); + } else if (ti.rfind("NEB=", 0) == 0) { + neb_ = parse_param(ti.substr(4), iObj.parameters, s.second); + } else if (ti.rfind("AREA=", 0) == 0) { + area_ = spread.spread_value( + parse_param(ti.substr(5), iObj.parameters, s.second), Spread::JJ, spr); + } else if (ti.rfind("IC=", 0) == 0) { + Ic_ = spread.spread_value( + parse_param(ti.substr(3), iObj.parameters, s.second), Spread::JJ, spr); + } else { + t.emplace_back(ti); } } + if (!temp_ && iObj.globalTemp) { + temp_ = iObj.globalTemp; + } + if (!neb_ && iObj.neB) { + neb_ = iObj.neB; + } // Set component timestep h_ = iObj.transSim.tstep(); at_ = iObj.argAnal; @@ -62,30 +82,8 @@ JJ::JJ( netlistInfo.label_ = s.first.at(0); // Add the label to the known labels list lm.emplace(s.first.at(0)); - tokens_t t; - // Identify the named parameters of the junction, leaving model as last token - for (auto i : s.first) { - if (i.find("AREA=") != std::string::npos) { - area_ = spread.spread_value( - parse_param(i.substr(i.find("AREA=") + 5), iObj.parameters, s.second), - Spread::JJ, spr); - } else if (i.find("IC=") != std::string::npos) { - Ic_ = spread.spread_value( - parse_param(i.substr(i.find("IC=") + 3), iObj.parameters, s.second), - Spread::JJ, spr); - } else if (i.find("SPREAD=") != std::string::npos) { - } else { - t.emplace_back(i); - } - } // Set the model for this JJ instance set_model(t, iObj.netlist.models_new, s.second); - // Get and set the phase offset from the model - pn1_ = model_.value().get_phaseOffset(); - // Set the transitional conductance value - gLarge_ = model_.value().get_criticalCurrent() / - (model_.value().get_criticalToNormalRatio() * - model_.value().get_deltaV()); // Set the phase constant if (at_ == AnalysisType::Voltage) { // If voltage mode set this to (3 * hbar) / (4 * h * eV) @@ -94,19 +92,6 @@ JJ::JJ( // If phase mode set this to (4 * h * eV) / (3 * hbar) phaseConst_ = (4 * h_ * Constants::EV) / (3 * Constants::HBAR); } - // Set the Del0 parameter - del0_ = 1.76 * Constants::BOLTZMANN * - model_.value().get_criticalTemperature(); - // Set the del parameter - del_ = del0_ * sqrt(cos((Constants::PI / 2) * - (model_.value().get_temperature() / - model_.value().get_criticalTemperature()) * - (model_.value().get_temperature() / - model_.value().get_criticalTemperature()))); - // Set the calculated Rn value - rncalc_ = ((Constants::PI * del_) / - (2 * Constants::EV * model_.value().get_criticalCurrent())) * - tanh(del_ / (2 * Constants::BOLTZMANN * model_.value().get_temperature())); // Set the node configuration type indexInfo.nodeConfig_ = ncon; // Set variable index and increment it @@ -117,24 +102,30 @@ JJ::JJ( set_node_indices(tokens_t(s.first.begin() + 1, s.first.begin() + 3), nm, nc); // Set the non zero, column index and row pointer vectors set_matrix_info(); + if (temp_) { + spAmp_ = Noise::determine_spectral_amplitude( + this->model_.r0(), temp_.value()); + Function tnoise; + tnoise.parse_function("NOISE(" + + Misc::precise_to_string(spAmp_.value()) + ", 0.0, " + + Misc::precise_to_string(1.0 / neb_.value()) + ")", iObj, s.second); + thermalNoise = tnoise; + } } double JJ::subgap_impedance() { // Set subgap impedance (1/R0) + (3C/2h) - return ((1 / model_.value().get_subgapResistance()) + - ((3.0 * model_.value().get_capacitance()) / (2.0 * h_))); + return ((1 / model_.r0()) + ((3.0 * model_.c()) / (2.0 * h_))); } double JJ::transient_impedance() { // Set transitional impedance (GL) + (3C/2h) - return (gLarge_ + - ((3.0 * model_.value().get_capacitance()) / (2.0 * h_))); + return (gLarge_ + ((3.0 * model_.c()) / (2.0 * h_))); } double JJ::normal_impedance() { // Set normal impedance (1/RN) + (3C/2h) - return ((1 / model_.value().get_normalResistance()) + - ((3.0 * model_.value().get_capacitance()) / (2.0 * h_))); + return ((1 / model_.rn()) + ((3.0 * model_.c()) / (2.0 * h_))); } void JJ::set_matrix_info() { @@ -185,61 +176,79 @@ void JJ::set_matrix_info() { void JJ::set_model( const tokens_t& t, const vector_pair_t& models, const string_o& subc) { + bool found = false; // Loop through all models for (auto& i : models) { // If the model name matches that of an identified model - if (i.first.get_modelName() == t.back()) { + if (i.first.modelName() == t.back()) { // If both models belong to a subcircuit and the subcircuit names match if ((i.second && subc) && (subc.value() == i.second.value())) { // Set the model to the exact identified model model_ = i.first; + found = true; break; - // JJ might be in a subcircuit but model in global scope + // JJ might be in a subcircuit but model in global scope6 } else if (!i.second) { // Set the model to the globally identified model model_ = i.first; + found = true; break; } } } // If no model was found - if (!model_) { + if (!found) { // Complain about it Errors::invalid_component_errors( ComponentErrors::MODEL_NOT_DEFINED, Misc::vector_to_string(t)); } + // Set the model critical current for this JJ instance + model_.ic(model_.ic() * area_); + if (model_.tDep()) { + // Set the Del0 parameter + del0_ = 1.76 * Constants::BOLTZMANN * model_.tc(); + // Set the del parameter + del_ = del0_ * sqrt(cos((Constants::PI / 2) * + (model_.t() / model_.tc()) * (model_.t() / model_.tc()))); + // Set the temperature dependent normal resistance + model_.rn(((Constants::PI * del_) / (2 * Constants::EV * model_.ic())) * + tanh(del_ / (2 * Constants::BOLTZMANN * model_.t()))); + } else { + // Set the model normal resistance for this JJ instance + model_.rn(model_.rn() / area_); + } // Change the area if ic was defined if (Ic_) { - area_ = Ic_.value() / model_.value().get_criticalCurrent(); + area_ = Ic_.value() / model_.ic(); } // Set the model capacitance for this JJ instance - model_.value().set_capacitance( - model_.value().get_capacitance() * area_); - // Set the model normal resistance for this JJ instance - model_.value().set_normalResistance( - model_.value().get_normalResistance() / area_); + model_.c(model_.c() * area_); // Set the model subgap resistance for this JJ instance - model_.value().set_subgapResistance( - model_.value().get_subgapResistance() / area_); - // Set the model critical current for this JJ instance - model_.value().set_criticalCurrent( - model_.value().get_criticalCurrent() * area_); + model_.r0(model_.r0() / area_); // Set the lower boundary for the transition region - lowerB_ = - model_.value().get_voltageGap() - 0.5 * model_.value().get_deltaV(); + lowerB_ = model_.vg() - 0.5 * model_.deltaV(); // Set the upper boundary for the transition region - upperB_ = - model_.value().get_voltageGap() + 0.5 * model_.value().get_deltaV(); + upperB_ = model_.vg() + 0.5 * model_.deltaV(); + // Set the transitional conductance value + gLarge_ = model_.ic() / (model_.icFct() * model_.deltaV()); + if (model_.rtype() == 0) { + model_.r0(model_.rn()); + } } // Update the value based on the matrix entry based on voltage value bool JJ::update_value(const double& v) { // Shorthand for the model - const Model& m = model_.value(); + const Model& m = model_; // If the absolute value of the voltage is less than lower bounds if (fabs(v) < lowerB_) { + // Set temperature resistance + if (this->temp_) { + thermalNoise.value().ampValues().at(0) = + Noise::determine_spectral_amplitude(this->model_.r0(), temp_.value()); + } // Set the transition current to 0 - transitionCurrent_ = 0.0; + it_ = 0.0; // If the non zero vector does not end with the subgap conductance if (matrixInfo.nonZeros_.back() != -1 / subgap_impedance()) { // Make it end with the subgap conductance @@ -255,9 +264,11 @@ bool JJ::update_value(const double& v) { // If the absolute value of the voltage is less than the upperbounds } else if (fabs(v) < upperB_) { // Set the transition current - transitionCurrent_ = lowerB_ * ((1 / m.get_subgapResistance()) - gLarge_); + it_ = lowerB_ * ((1 / m.r0()) - gLarge_); // If the voltage is negative, current must be negative - if (v < 0) transitionCurrent_ = -transitionCurrent_; + if (v < 0) { + it_ = -it_; + } // If the back of the non zero vector is not the transition conductance if (matrixInfo.nonZeros_.back() != -1 / transient_impedance()) { // Set it to the transition conductance @@ -272,8 +283,13 @@ bool JJ::update_value(const double& v) { } // If neither of the above then it must be in normal operating region } else { + // Set temperature resistance + if (this->temp_) { + thermalNoise.value().ampValues().at(0) = + Noise::determine_spectral_amplitude(this->model_.rn(), temp_.value()); + } // Reset the transition current, transition has passed. - transitionCurrent_ = 0.0; + it_ = 0.0; // If the back of the non zero vector is not the normal conductance if (matrixInfo.nonZeros_.back() != -1 / normal_impedance()) { // Set it to the normal conductance @@ -289,23 +305,4 @@ bool JJ::update_value(const double& v) { } // This should never happen, only here to appease the return type return false; -} - -// Update timestep based on a scalar factor i.e 0.5 for half the timestep -void JJ::update_timestep(const double& factor) { - h_ = h_ * factor; - if (state_ == 0) { - matrixInfo.nonZeros_.back() = -1 / subgap_impedance(); - } else if (state_ == 1) { - matrixInfo.nonZeros_.back() = -1 / transient_impedance(); - } else if (state_ == 2) { - matrixInfo.nonZeros_.back() = -1 / normal_impedance(); - } - if (at_ == AnalysisType::Voltage) { - matrixInfo.nonZeros_.at(hDepPos_) = - (1 / factor) * matrixInfo.nonZeros_.at(hDepPos_); - } else if (at_ == AnalysisType::Phase) { - matrixInfo.nonZeros_.at(hDepPos_) = - factor * matrixInfo.nonZeros_.at(hDepPos_); - } } \ No newline at end of file diff --git a/src/LUSolve.cpp b/src/LUSolve.cpp index b13ce6d..05c0542 100644 --- a/src/LUSolve.cpp +++ b/src/LUSolve.cpp @@ -63,9 +63,9 @@ void LUSolve::solve(std::vector& x) { if (!constructed) { ABORT("Preconditioner not constructed."); } - int x_size = x.size(); - DNformat LHSstore = { x_size, &x.front() }; - SuperMatrix LHSmat = { SLU_DN, SLU_D, SLU_GE, x_size, 1, &LHSstore }; + DNformat LHSstore = { static_cast(x.size()), &x.front() }; + SuperMatrix LHSmat = { SLU_DN, SLU_D, SLU_GE, static_cast(x.size()), + 1, &LHSstore }; dgstrs(trans, &L, &U, perm_c, perm_r, &LHSmat, &stat, &info); } diff --git a/src/Matrix.cpp b/src/Matrix.cpp index 5ca912a..d6d6fa0 100644 --- a/src/Matrix.cpp +++ b/src/Matrix.cpp @@ -16,16 +16,33 @@ using namespace JoSIM; void Matrix::create_matrix(Input& iObj) { + while (needsTR_) { + // Do matrix setup + setup(iObj); + // Create the components + create_components(iObj); + if (needsTR_) { + reduce_step(iObj); + } + } + // Handle mutual inductances + handle_mutual_inductance(iObj); + // Create the compressed storage row format required for simulation + create_csr(); +} + +void Matrix::setup(Input& iObj) { + needsTR_ = false; spread.get_spreads(iObj); - // Create a seperate thread that will be used for printing creation progress - std::thread printingThread; + Noise::determine_global_temperature(iObj); + Noise::determine_noise_effective_bandwidth(iObj); // Create a node counter variable int nodeCounter = 0; // Variables to store node configs since they are already identified here - std::vector nodeConfig( + nodeConfig.resize( iObj.netlist.expNetlist.size(), NodeConfig::GND); - // Optional, only available when 4 node components are detected - std::vector> nodeConfig2( + // Only available when 4 node components are detected + nodeConfig2.resize( iObj.netlist.expNetlist.size(), NodeConfig::GND); int cc = 0; for (auto& i : iObj.netlist.expNetlist) { @@ -40,13 +57,13 @@ void Matrix::create_matrix(Input& iObj) { Misc::vector_to_string(i.first)); } // Create a node map that maps the node names to numbers - if (i.first.at(1).find("GND") == + if (i.first.at(1).find("GND") == std::string::npos && i.first.at(1) != "0") { // Add the first node to the map if not ground if (nm.count(i.first.at(1)) == 0) nm[i.first.at(1)] = nodeCounter++; nodeConfig.at(cc) = NodeConfig::POSGND; } - if (i.first.at(2).find("GND") == + if (i.first.at(2).find("GND") == std::string::npos && i.first.at(2) != "0") { // Add the second node to the map if not ground if (nm.count(i.first.at(2)) == 0) nm[i.first.at(2)] = nodeCounter++; @@ -57,7 +74,7 @@ void Matrix::create_matrix(Input& iObj) { } } // If the device is as 4 node device - if (std::string("EFGHT").find(i.first.front().at(0)) != + if (std::string("EFGHT").find(i.first.front().at(0)) != std::string::npos) { // Ensure the device has at least 6 parts if (i.first.size() < 6) { @@ -90,6 +107,10 @@ void Matrix::create_matrix(Input& iObj) { nc.resize(nm.size()); // Set the index to the first branch current to the size of the node map branchIndex = nm.size(); +} + +void Matrix::create_components(Input& iObj) { + int cc = 0; ProgressBar bar; if (!iObj.argMin) { bar.create_thread(); @@ -99,8 +120,6 @@ void Matrix::create_matrix(Input& iObj) { bar.set_status_text("Creating Matrix"); bar.set_total((float)iObj.netlist.expNetlist.size()); } - // Counter for progress report - cc = 0; // Loop through all the components in the netlist for (auto it = iObj.netlist.expNetlist.begin(); it != iObj.netlist.expNetlist.end(); it++) { @@ -198,6 +217,11 @@ void Matrix::create_matrix(Input& iObj) { i, nodeConfig.at(cc), nodeConfig2.at(cc), nm, lm, nc, iObj.parameters, iObj.argAnal, iObj.transSim.tstep(), branchIndex)); + if (std::get(components.devices.back()).timestepDelay_ + <= 4) { + needsTR_ = true; + return; + } // Store the transmission line component list index for reference components.txIndices.emplace_back(components.devices.size() - 1); break; @@ -255,6 +279,10 @@ void Matrix::create_matrix(Input& iObj) { bar.complete(); std::cout << "\n"; } +} + +void Matrix::handle_mutual_inductance(Input& iObj) { + int cc = 0; ProgressBar bar2; if (!iObj.argMin && components.mutualinductances.size() != 0) { bar2.create_thread(); @@ -264,8 +292,6 @@ void Matrix::create_matrix(Input& iObj) { bar2.set_status_text("Adding Mutual Inductances"); bar2.set_total((float)iObj.netlist.expNetlist.size()); } - // Counter for progress report - cc = 0; // Loop through all identified mutual inductances for (const auto& s : components.mutualinductances) { // If not minimal printing @@ -329,14 +355,29 @@ void Matrix::create_matrix(Input& iObj) { } ++cc; } - // Create the compressed storage row format required for simulation - create_csr(); // Led the user know the matrix creation is complete if (!iObj.argMin && components.mutualinductances.size() != 0) { bar2.complete(); std::cout << "\n"; } +} +void Matrix::reduce_step(Input& iObj) { + iObj.transSim.tstep(iObj.transSim.tstep() / 2); + nodeConfig.clear(); + nodeConfig2.clear(); + Components newComponents; + components = newComponents; + nm.clear(); + nodeconnections newNC; + nc = newNC; + lm.clear(); + branchIndex = 0; + nz.clear(); + ci.clear(); + rp.clear(); + relevantIndices.clear(); + relevantTraces.clear(); } void Matrix::create_csr() { diff --git a/src/Misc.cpp b/src/Misc.cpp index 5a44f6b..3416569 100644 --- a/src/Misc.cpp +++ b/src/Misc.cpp @@ -96,9 +96,11 @@ tokens_t Misc::tokenize( } // If trim spaces is enabled if (trimSpaces) { - // Remove trailing, leading and duplicate spaces between tokens - tokens.back() = - std::regex_replace(tokens.back(), std::regex("^ +| +$|( ) +"), "$1"); + if (!tokens.empty()) { + // Remove trailing, leading and duplicate spaces between tokens + tokens.back() = + std::regex_replace(tokens.back(), std::regex("^ +| +$|( ) +"), "$1"); + } } lastPos = pos + 1; ++counter; diff --git a/src/Model.cpp b/src/Model.cpp index 5869862..a0c8c40 100644 --- a/src/Model.cpp +++ b/src/Model.cpp @@ -2,8 +2,8 @@ // This code is licensed under MIT license (see LICENSE for details) #include "JoSIM/Model.hpp" -#include "JoSIM/Misc.hpp" #include "JoSIM/Errors.hpp" +#include "JoSIM/Misc.hpp" #include "JoSIM/TypeDefines.hpp" using namespace JoSIM; @@ -19,7 +19,7 @@ void Model::parse_model( // Create a model that will be stored Model temp; // The second token is the model label - temp.set_modelName(s.first.at(1)); + temp.modelName(s.first.at(1)); // The third token needs to start with "JJ" for this to be valid if (s.first.at(2).compare(0, 2, "JJ") != 0) { Errors::model_errors( @@ -34,39 +34,49 @@ void Model::parse_model( Errors::model_errors( ModelErrors::BAD_MODEL_DEFINITION, Misc::vector_to_string(s.first)); } + double value = 0.0; // Loop through the parameter tokens for (int i = 0; i < tokens.size(); i += 2) { // Every even odd token should be a value (otherwise complain) - double value = parse_param(tokens.at(i + 1), p, s.second); + value = parse_param(tokens.at(i + 1), p, s.second); if (std::isnan(value)) { Errors::model_errors( ModelErrors::BAD_MODEL_DEFINITION, Misc::vector_to_string(s.first)); } // Assign the relevant model parameters if (tokens.at(i) == "VG" || tokens.at(i) == "VGAP") { - temp.set_voltageGap(value); + temp.vg(value); } else if (tokens.at(i) == "IC" || tokens.at(i) == "ICRIT") { - temp.set_criticalCurrent(value); + temp.ic(value); } else if (tokens.at(i) == "RTYPE") { - temp.set_resistanceType((int)value); + temp.rtype((int)value); } else if (tokens.at(i) == "RN") { - temp.set_normalResistance(value); + temp.rn(value); } else if (tokens.at(i) == "R0") { - temp.set_subgapResistance(value); + temp.r0(value); } else if (tokens.at(i) == "CAP" || tokens.at(i) == "C") { - temp.set_capacitance(value); + temp.c(value); } else if (tokens.at(i) == "T") { - temp.set_temperature(value); + temp.t(value); + temp.tDep(true); } else if (tokens.at(i) == "TC") { - temp.set_criticalTemperature(value); + temp.tc(value); + temp.tDep(true); } else if (tokens.at(i) == "DELV") { - temp.set_deltaV(value); + temp.deltaV(value); } else if (tokens.at(i) == "D") { - temp.set_transparency(value); + temp.d(value); + temp.tDep(true); } else if (tokens.at(i) == "ICFACT" || tokens.at(i) == "ICFCT") { - temp.set_criticalToNormalRatio(value); + temp.icFct(value); } else if (tokens.at(i) == "PHI") { - temp.set_phaseOffset(value); + temp.phiOff(value); + } else { + // Incompatible parameter + Errors::model_errors(ModelErrors::PARAM_TYPE_ERROR, + Misc::vector_to_string( + tokens_t{ Misc::vector_to_string(s.first), + "\nThe parameter: ", tokens.at(i) })); } } diff --git a/src/Netlist.cpp b/src/Netlist.cpp index c8626c7..66241b5 100644 --- a/src/Netlist.cpp +++ b/src/Netlist.cpp @@ -13,28 +13,30 @@ using namespace JoSIM; void Netlist::id_io_subc_label( - const tokens_t& lineTokens, tokens_t& io, + const tokens_t& lineTokens, tokens_t& io, s_map& params, std::string& subcktName, std::string& label, const std::unordered_map& subcircuits) { // Id the label label = lineTokens.front(); + // Found flag + bool found = false; // Check the convention // At this point in the program all subcircuits should have been identified // Thus we can determine the convention from code - // If the token right after the label exists in the subcircuits map - if (subcircuits.count(lineTokens.at(1)) != 0) { - // This is the subcircuit name - subcktName = lineTokens.at(1); - // Assign the IO - io.assign(lineTokens.begin() + 2, lineTokens.end()); - // Else if the last token exists in the subcircuit map - } else if (subcircuits.count(lineTokens.back()) != 0) { - // Then this is the subcircuit name - subcktName = lineTokens.back(); - // Assign the IO - io.assign(lineTokens.begin() + 1, lineTokens.end() - 1); - // Else if the neither then this subcircuit surely does not exist - } else { + for (int i = 1; i < lineTokens.size(); ++i) { + if (subcircuits.count(lineTokens.at(i)) != 0) { + found = true; + subcktName = lineTokens.at(i); + } else if (lineTokens.at(i).find('=') != std::string::npos) { + std::string param, value; + param = lineTokens.at(i).substr(0, lineTokens.at(i).find('=')); + value = lineTokens.at(i).substr(lineTokens.at(i).find('=') + 1); + params[param] = value; + } else { + io.emplace_back(lineTokens.at(i)); + } + } + if (!found) { Errors::input_errors( InputErrors::UNKNOWN_SUBCKT, Misc::vector_to_string(lineTokens)); } @@ -59,12 +61,20 @@ bool Netlist::rename_io_nodes( return false; } -void Netlist::expand_io( - Subcircuit& subc, const tokens_t& io, const std::string& label) { +void Netlist::expand_io(Subcircuit& subc, s_map& params, const tokens_t& io, + const std::string& label) { + // Sanity check + if (io.size() != subc.io.size()) { + Errors::input_errors(InputErrors::IO_MISMATCH, label); + } // Loop through the identified subcircuit for (int k = 0; k < subc.lines.size(); ++k) { // Set shorthand for long variable name tokens_t& tokens = subc.lines.at(k).first; + // Check for value parameterization + if (params.count(tokens.front()) != 0) { + insert_parameter(tokens, params); + } // Append the label of the parent to the subcircuit label tokens.at(0) = tokens.at(0) + "|" + label; // Determine amount of nodes to process @@ -87,9 +97,33 @@ void Netlist::expand_io( } } +void Netlist::insert_parameter(tokens_t& t, s_map& params) { + std::string twonode = "RCLK"; + std::string fournode = "EFHGT"; + if (twonode.find(t.front().front()) != std::string::npos) { + t.at(3) = params.at(t.front()); + } else if (fournode.find(t.front().front()) != std::string::npos) { + if (t.front().front() == 'T') { + for (auto& i : t) { + if (i.find("TD=") != std::string::npos) { + i = "TD=" + params.at(t.front()); + } + } + } else { + t.at(5) = params.at(t.front()); + } + } else if (t.front().front() == 'B') { + for (auto& i : t) { + if (i.find("AREA=") != std::string::npos) { + i = "AREA=" + params.at(t.front()); + } else if (i.find("IC=") != std::string::npos) { + i = "IC=" + params.at(t.front()); + } + } + } +} + void Netlist::expand_subcircuits() { - // Variable to store io - tokens_t io; // std::vector> moddedLines; std::string subcktName, label; // Loop through subcircuits, line by line @@ -133,14 +167,17 @@ void Netlist::expand_subcircuits() { tokens_t& subcCurrentLine = subcircuit.lines.at(j).first; // If the line denotes a subcircuit if (subcCurrentLine.front().at(0) == 'X') { + // Variable to store io and parameters + tokens_t io; + s_map params; id_io_subc_label( - subcCurrentLine, io, subcktName, label, subcircuits); + subcCurrentLine, io, params, subcktName, label, subcircuits); // Create a copy of the subircuit for this instance Subcircuit subc = subcircuits.at(subcktName); // If not nested if (!subc.containsSubckt) { // Expand the IO nodes of the subcircuit for this instance - expand_io(subc, io, label); + expand_io(subc, params, io, label); // Erase the current line (tokens) subcircuit.lines.erase(subcircuit.lines.begin() + j); // Insert the subcircuit token lines at the erased line position @@ -168,8 +205,6 @@ void Netlist::expand_subcircuits() { } void Netlist::expand_maindesign() { - // Variable to store io - tokens_t io; // std::vector> moddedLines; ProgressBar bar; if (!argMin) { @@ -190,11 +225,15 @@ void Netlist::expand_maindesign() { } // If the line denotes a subcircuit if (maindesign.at(i).front().at(0) == 'X') { - id_io_subc_label(maindesign.at(i), io, subcktName, label, subcircuits); + // Variable to store io and parameters + tokens_t io; + s_map params; + id_io_subc_label(maindesign.at(i), io, params, + subcktName, label, subcircuits); // Copy of subcircuit for this instance Subcircuit subc = subcircuits.at(subcktName); // Expand the appropriate IO nodes of the subcircuit for this instance - expand_io(subc, io, label); + expand_io(subc, params, io, label); // Add the expanded subcircuit lines to the expanded netlist expNetlist.insert(expNetlist.end(), subc.lines.begin(), subc.lines.end()); diff --git a/src/Noise.cpp b/src/Noise.cpp index d3b8ef2..7732c59 100644 --- a/src/Noise.cpp +++ b/src/Noise.cpp @@ -27,79 +27,7 @@ void JoSIM::Noise::determine_noise_effective_bandwidth(Input& iObj) { double JoSIM::Noise::determine_spectral_amplitude(const double& R, const double& T) { - // spectral amplitute = sqrt(4 * kB * T * B / R) + // spectral amplitute = sqrt(4 * kB * T / R) double spAmp = sqrt((4 * Constants::BOLTZMANN * T) / R); return spAmp; } - -std::pair -JoSIM::Noise::create_resistive_current_noise(Input& iObj, const double& R, - const double& T, const double& B, - std::vector>::iterator& i) { - std::pair noiseSource; - noiseSource.first.emplace_back("INOISE_" + i->first.front()); - noiseSource.first.emplace_back(i->first.at(1)); - noiseSource.first.emplace_back(i->first.at(2)); - noiseSource.first.emplace_back("NOISE(" + Misc::precise_to_string( - determine_spectral_amplitude(R, T))); - noiseSource.first.emplace_back("0"); - noiseSource.first.emplace_back( - Misc::precise_to_string(1 / B) + ")"); - noiseSource.second = i->second; - return noiseSource; -} - -void JoSIM::Noise::add_noise_sources(Input& iObj) { - determine_global_temperature(iObj); - determine_noise_effective_bandwidth(iObj); - std::vector> noises; - for (auto it = iObj.netlist.expNetlist.begin(); - it != iObj.netlist.expNetlist.end(); ++it) { - if (it->first.front().at(0) == 'R') { - std::optional T; - double B = 0, R = 0; - if (!iObj.neB) { - B = 1E12; - } else { - B = iObj.neB.value(); - } - if (iObj.globalTemp) { - T = iObj.globalTemp.value(); - } - if (it->first.size() > 4) { - tokens_t t(it->first.begin() + 4, it->first.end()); - std::string temp = Misc::vector_to_string(t, " "); - t = Misc::tokenize(temp, " ="); - for (auto i = t.begin(); i < t.end(); ++i) { - if (*i == "TEMP") { - if (i + 1 == t.end() || *(i + 1) == "NEB") { - Errors::invalid_component_errors(ComponentErrors::RES_ERROR, - Misc::vector_to_string(it->first)); - } else { - T = parse_param(*(i + 1), iObj.parameters, it->second); - } - } - if (*i == "NEB") { - if (i + 1 == t.end() || *(i + 1) == "TEMP") { - Errors::invalid_component_errors(ComponentErrors::RES_ERROR, - Misc::vector_to_string(it->first)); - } else { - B = parse_param(*(i + 1), iObj.parameters, it->second); - } - } - } - R = parse_param(*(it->first.begin() + 3), iObj.parameters, it->second); - } else { - R = parse_param(it->first.back(), iObj.parameters, it->second); - } - if (T) { - noises.emplace_back( - create_resistive_current_noise(iObj, R, T.value(), B, it)); - } - } - } - if (noises.size() != 0) { - iObj.netlist.expNetlist.insert(iObj.netlist.expNetlist.end(), - noises.begin(), noises.end()); - } -} diff --git a/src/Parameters.cpp b/src/Parameters.cpp index dc156a4..775e602 100644 --- a/src/Parameters.cpp +++ b/src/Parameters.cpp @@ -73,6 +73,10 @@ double JoSIM::parse_param( bool single) { // Initialize the expression to evaluate std::string expToEval = expr; + // Sanity check, prepend 0 to a string where the value is eg. .5 to vorm 0.5 + if (expToEval.front() == '.') { + expToEval = "0" + expToEval; + } // Remove any and all whitespace characters expToEval.erase(std::remove_if(expToEval.begin(), expToEval.end(), isspace), expToEval.end()); @@ -476,41 +480,33 @@ void JoSIM::expand_inline_parameters( Errors::parsing_errors( ParsingErrors::INVALID_DECLARATION, Misc::vector_to_string(s.first)); } - // Create a temporary vector to split the token on '{' - tokens_t temp = - Misc::tokenize(s.first.at(oPos.value()), "{", true, true, 1); - // Erase the original token containing the '{' - s.first.erase(s.first.begin() + oPos.value()); - // Insert at the previous token position the temporary vector - s.first.insert(s.first.begin() + oPos.value(), temp.begin(), temp.end()); - // If the emporary token is of size 2. e.g. "td={4" becomes [td=][4] - if (temp.size() == 2) { - // Increment the position to the begining of the expression - oPos = oPos.value() + 1; + if (cPos.value() != oPos.value()) { + s.first.at(oPos.value()) = Misc::vector_to_string( + tokens_t{ s.first.begin() + oPos.value(), + s.first.begin() + cPos.value() + 1 }, ""); + s.first.erase(s.first.begin() + oPos.value() + 1, + s.first.begin() + cPos.value() + 1); } - // Split the closing parenthesis into tokens on '}' - temp = - Misc::tokenize(s.first.at(cPos.value()), "}", true, true, 1); - // Erase the original closing token - s.first.erase(s.first.begin() + cPos.value()); - // Insert the temporary vector at the original position - s.first.insert(s.first.begin() + cPos.value(), temp.begin(), temp.end()); - // Parse the identified expression - double value = parse_param( - Misc::vector_to_string(tokens_t(s.first.begin() + oPos.value(), - s.first.begin() + cPos.value()), ""), parameters, s.second, false); + tokens_t temp = Misc::tokenize(s.first.at(oPos.value()), "{}", true, true, 2); + double value = parse_param(temp.back(), parameters, s.second, false); // If the returned value is not NaN - if (!std::isnan(value)) { + if (std::isnan(value)) { // Complain of invalid parameter expression Errors::parsing_errors( ParsingErrors::UNIDENTIFIED_PART, Misc::vector_to_string(s.first)); } // Erase the expression part of the tokens s.first.erase( - s.first.begin() + oPos.value(), s.first.begin() + cPos.value()); - // Insert the double value in the place of the expression - s.first.insert( - s.first.begin() + oPos.value(), Misc::precise_to_string(value)); + s.first.begin() + oPos.value(), s.first.begin() + oPos.value() + 1); + if (temp.size() > 1) { + // Insert the double value in the place of the expression + s.first.insert(s.first.begin() + oPos.value(), + Misc::vector_to_string( + tokens_t{ temp.front(), Misc::precise_to_string(value) }, "")); + } else { + s.first.insert(s.first.begin() + oPos.value(), + Misc::precise_to_string(value)); + } } // If an expression closing was found if (cPos) { diff --git a/src/Resistor.cpp b/src/Resistor.cpp index 1e5d214..017a9c4 100644 --- a/src/Resistor.cpp +++ b/src/Resistor.cpp @@ -5,6 +5,7 @@ #include "JoSIM/Misc.hpp" #include "JoSIM/Errors.hpp" #include "JoSIM/Constants.hpp" +#include "JoSIM/Noise.hpp" #include @@ -31,12 +32,22 @@ Resistor::Resistor( const nodemap& nm, std::unordered_set& lm, nodeconnections& nc, Input& iObj, Spread& spread, int& bi) { double spr = 1.0; - for (auto i : s.first) { - if (i.find("SPREAD=") != std::string::npos) { - spr = - parse_param(i.substr(i.find("SPREAD=") + 7), iObj.parameters, s.second); + for (int i = 4; i < s.first.size(); ++i) { + auto& t = s.first.at(i); + if (t.rfind("SPREAD=", 0) == 0) { + spr = parse_param(t.substr(7), iObj.parameters, s.second); + } else if (t.rfind("TEMP=", 0) == 0) { + temp_ = parse_param(t.substr(5), iObj.parameters, s.second); + } else if (t.rfind("NEB=", 0) == 0) { + neb_ = parse_param(t.substr(4), iObj.parameters, s.second); } } + if (!temp_ && iObj.globalTemp) { + temp_ = iObj.globalTemp; + } + if (!neb_ && iObj.neB) { + neb_ = iObj.neB; + } at_ = iObj.argAnal; // Check if the label has already been defined if (lm.count(s.first.at(0)) != 0) { @@ -70,6 +81,15 @@ Resistor::Resistor( -((2.0 * iObj.transSim.tstep()) / 3.0) * (netlistInfo.value_ / Constants::SIGMA)); } + if (temp_) { + spAmp_ = Noise::determine_spectral_amplitude( + netlistInfo.value_, temp_.value()); + Function tnoise; + tnoise.parse_function("NOISE(" + + Misc::precise_to_string(spAmp_.value()) + ", 0.0, " + + Misc::precise_to_string(1.0 / neb_.value()) + ")", iObj, s.second); + thermalNoise = tnoise; + } } // Update timestep based on a scalar factor i.e 0.5 for half the timestep diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 712f8b0..f715bae 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -15,6 +15,45 @@ using namespace JoSIM; Simulation::Simulation(Input &iObj, Matrix &mObj) { + while (needsTR_) { + // Do generic simulation setup for given step size + setup(iObj, mObj); + // Do solver setup + if (SLU) { + // SLU setup + lu.create_matrix(mObj.rp.size() - 1, mObj.nz, mObj.ci, mObj.rp); + lu.factorize(); + } else { + // KLU setup + simOK_ = klu_l_defaults(&Common_); + assert(simOK_); + Symbolic_ = klu_l_analyze( + mObj.rp.size() - 1, &mObj.rp.front(), &mObj.ci.front(), &Common_); + Numeric_ = klu_l_factor( + &mObj.rp.front(), &mObj.ci.front(), &mObj.nz.front(), + Symbolic_, &Common_); + } + + // Run transient simulation + trans_sim(mObj); + // If step size is too large, reduce and try again + if (needsTR_) { + reduce_step(iObj, mObj); + } + + // Do solver cleanup + if (SLU) { + // SLU cleanup + lu.free(); + } else { + // KLU cleanup + klu_l_free_symbolic(&Symbolic_, &Common_); + klu_l_free_numeric(&Numeric_, &Common_); + } + } +} + +void Simulation::setup(Input& iObj, Matrix& mObj) { // Simulation setup SLU = iObj.SLU; simSize_ = iObj.transSim.simsize(); @@ -25,38 +64,17 @@ Simulation::Simulation(Input &iObj, Matrix &mObj) { stepSize_ = iObj.transSim.tstep(); prstep_ = iObj.transSim.prstep(); prstart_ = iObj.transSim.prstart(); + startup_ = iObj.transSim.startup(); + x_.clear(); x_.resize(mObj.branchIndex, 0.0); - if(!mObj.relevantTraces.empty()) { + if (!mObj.relevantTraces.empty()) { results.xVector.resize(mObj.branchIndex); - for (const auto &i : mObj.relevantIndices) { + for (const auto& i : mObj.relevantIndices) { results.xVector.at(i).emplace(); } } else { results.xVector.resize(mObj.branchIndex, std::vector(0)); } - if (SLU) { - // SLU setup - lu.create_matrix(mObj.rp.size() - 1, mObj.nz, mObj.ci, mObj.rp); - lu.factorize(); - // Run transient simulation - trans_sim(mObj); - // SLU cleanup - lu.free(); - } else { - // KLU setup - simOK_ = klu_l_defaults(&Common_); - assert(simOK_); - Symbolic_ = klu_l_analyze( - mObj.rp.size() - 1, &mObj.rp.front(), &mObj.ci.front(), &Common_); - Numeric_ = klu_l_factor( - &mObj.rp.front(), &mObj.ci.front(), &mObj.nz.front(), - Symbolic_, &Common_); - // Run transient simulation - trans_sim(mObj); - // KLU cleanup - klu_l_free_symbolic(&Symbolic_, &Common_); - klu_l_free_numeric(&Numeric_, &Common_); - } } void Simulation::trans_sim(Matrix &mObj) { @@ -73,6 +91,28 @@ void Simulation::trans_sim(Matrix &mObj) { } // Initialize the b matrix b_.resize(mObj.rp.size(), 0.0); + if (startup_) { + // Stabilize the simulation before starting at t=0 + int startup = 2 * pow(10, (abs(log10(stepSize_)) - 12) * 2 + 1); + for (int i = -startup; i < 0; ++i) { + double step = i * stepSize_; + // Setup the b matrix + setup_b(mObj, i, i * stepSize_); + if (needsTR_) return; + // Assign x_prev the new b + x_ = b_; + // Solve Ax=b, storing the results in x_ + if (SLU) { + lu.solve(x_); + } else { + simOK_ = klu_l_tsolve( + Symbolic_, Numeric_, mObj.rp.size() - 1, 1, &x_.front(), &Common_); + // If anything is a amiss, complain about it + if (!simOK_) Errors::simulation_errors( + SimulationErrors::MATRIX_SINGULAR); + } + } + } // Start the simulation loop for(int i = 0; i < simSize_; ++i) { double step = i * stepSize_; @@ -82,6 +122,7 @@ void Simulation::trans_sim(Matrix &mObj) { } // Setup the b matrix setup_b(mObj, i, i * stepSize_); + if (needsTR_) return; // Assign x_prev the new b x_ = b_; // Solve Ax=b, storing the results in x_ @@ -109,6 +150,26 @@ void Simulation::trans_sim(Matrix &mObj) { } } +void Simulation::reduce_step(Input& iObj, Matrix& mObj) { + iObj.transSim.tstep(iObj.transSim.tstep() / 2); + bool tempMinOut = iObj.argMin; + if (!iObj.argMin) iObj.argMin = true; + Matrix newmObj; + mObj = newmObj; + // Create the matrix in csr format + for (auto& i : mObj.sourcegen) { + i.clearMisc(); + } + mObj.create_matrix(iObj); + find_relevant_traces(iObj, mObj); + //// Dump expanded Netlist since it is no longer needed + //iObj.netlist.expNetlist.clear(); + //iObj.netlist.expNetlist.shrink_to_fit(); + if (!tempMinOut) iObj.argMin = tempMinOut; + results.xVector.clear(); + results.timeAxis.clear(); +} + void Simulation::setup_b( Matrix &mObj, int i, double step, double factor) { // Clear b matrix and reset @@ -133,7 +194,7 @@ void Simulation::setup_b( // Handle current sources handle_cs(mObj, step, i); // Handle resistors - handle_resistors(mObj); + handle_resistors(mObj, step); // Handle inductors handle_inductors(mObj, factor); // Handle capacitors @@ -150,83 +211,6 @@ void Simulation::setup_b( handle_tx(mObj, i, step); } -void Simulation::reduce_step( - Matrix &mObj, double factor, - int &stepCount, double currentStep) { - // Backup the current nonzeros - mObj.nz_orig = mObj.nz; - // Restore the previous x - x_ = xn3_; - // Remove the stored results for the larger timesteps - for(auto &i : results.xVector){ - if(i) { - i.value().pop_back(); - i.value().pop_back(); - } - } - results.timeAxis.pop_back(); - results.timeAxis.pop_back(); - // Determine how many steps to take using the reduced time step - int smallSteps = static_cast(20E-12 / (factor * stepSize_)); - // Update the non-zeros of each component to reflect the smaller timestep - for (auto &j : mObj.components.devices) { - BasicComponent &x = std::visit( - [](auto& x) -> BasicComponent &{ return x; }, - j); - x.update_timestep(factor); - x.step_back(); - } - // Recreate the non-zero matrix for the simulation - mObj.create_nz(); - // Do a new LU decomposition - klu_l_free_numeric(&Numeric_, &Common_); - Numeric_ = klu_l_factor( - &mObj.rp.front(), &mObj.ci.front(), &mObj.nz.front(), - Symbolic_, &Common_); - for(int i = 0; i < smallSteps; ++i) { - double smallStep = currentStep + i * (factor * stepSize_); - // Setup the b matrix - setup_b( - mObj, stepCount, smallStep, factor); - // Assign x_prev the new b - x_ = b_; - // Solve Ax=b, storing the results in x_ - simOK_ = klu_l_tsolve( - Symbolic_, Numeric_, mObj.rp.size() - 1, 1, &x_.front(), &Common_); - // If anything is a amiss, complain about it - if (!simOK_) Errors::simulation_errors(SimulationErrors::MATRIX_SINGULAR); - float test = smallStep / stepSize_; - // if(test == stepCount) { - // Store results (only requested, to prevent massive memory usage) - for(int j = 0; j < results.xVector.size(); ++j) { - if(smallStep >= prstart_) { - if(results.xVector.at(j)) { - results.xVector.at(j).value().emplace_back(x_.at(j)); - } - } - } - // Store the time step - if(smallStep >= prstart_) results.timeAxis.emplace_back(smallStep); - if(test == stepCount) { - stepCount++; - } - } - // Restore the timestep of all components to the correct values - for (auto &j : mObj.components.devices) { - BasicComponent &x = std::visit( - [](auto& x) -> BasicComponent &{ return x; }, - j); - x.update_timestep(1/factor); - } - // Recreate the non-zero matrix for the simulation - mObj.create_nz(); - // Do a new LU decomposition - klu_l_free_numeric(&Numeric_, &Common_); - Numeric_ = klu_l_factor( - &mObj.rp.front(), &mObj.ci.front(), &mObj.nz.front(), - Symbolic_, &Common_); -} - void Simulation::handle_cs(Matrix &mObj, double &step, const int &i) { for (const auto &j : mObj.components.currentsources) { @@ -245,15 +229,29 @@ void Simulation::handle_cs(Matrix &mObj, double &step, const int &i) { } } -void Simulation::handle_resistors(Matrix &mObj) { +void Simulation::handle_resistors(Matrix &mObj, double &step) { for (const auto &j : mObj.components.resistorIndices) { auto &temp = std::get(mObj.components.devices.at(j)); NodeConfig &nc = temp.indexInfo.nodeConfig_; if(nc == NodeConfig::POSGND) { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.posIndex_.value()) -= + temp.thermalNoise.value().value(step); + } temp.pn1_ = (x_.at(temp.indexInfo.posIndex_.value())); } else if(nc == NodeConfig::GNDNEG) { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.negIndex_.value()) += + temp.thermalNoise.value().value(step); + } temp.pn1_ = (-x_.at(temp.indexInfo.negIndex_.value())); } else { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.posIndex_.value()) -= + temp.thermalNoise.value().value(step); + b_.at(temp.indexInfo.negIndex_.value()) += + temp.thermalNoise.value().value(step); + } temp.pn1_ = (x_.at(temp.indexInfo.posIndex_.value()) - x_.at(temp.indexInfo.negIndex_.value())); } @@ -330,10 +328,24 @@ void Simulation::handle_jj( auto &temp = std::get(mObj.components.devices.at(j)); const auto &model = temp.model_; if(temp.indexInfo.posIndex_ && !temp.indexInfo.negIndex_) { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.posIndex_.value()) -= + temp.thermalNoise.value().value(step); + } temp.pn1_ = (x_.at(temp.indexInfo.posIndex_.value())); } else if(!temp.indexInfo.posIndex_ && temp.indexInfo.negIndex_) { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.negIndex_.value()) += + temp.thermalNoise.value().value(step); + } temp.pn1_ = (-x_.at(temp.indexInfo.negIndex_.value())); } else { + if (temp.thermalNoise) { + b_.at(temp.indexInfo.posIndex_.value()) -= + temp.thermalNoise.value().value(step); + b_.at(temp.indexInfo.negIndex_.value()) += + temp.thermalNoise.value().value(step); + } temp.pn1_ = (x_.at(temp.indexInfo.posIndex_.value()) - x_.at(temp.indexInfo.negIndex_.value())); } @@ -352,21 +364,21 @@ void Simulation::handle_jj( // Phase guess (P0) temp.phi0_ = (4.0/3.0) * temp.pn1_ - (1.0/3.0) * temp.pn2_ + ((1.0 / Constants::SIGMA) * - ((2.0 * (stepSize_ * factor)) / 3.0)) * v0; + ((2.0 * (stepSize_)) / 3.0)) * v0; // Ensure timestep is not too large if ((double)i/(double)simSize_ > 0.01) { - if (abs(temp.phi0_ - temp.pn1_) > (0.25 * 2 * Constants::PI)) { - // needsTR_ = true; - // return; - Errors::simulation_errors( - SimulationErrors::PHASEGUESS_TOO_LARGE, temp.netlistInfo.label_); + if (abs(temp.phi0_ - temp.pn1_) > (0.20 * 2 * Constants::PI)) { + needsTR_ = true; + return; + /*Errors::simulation_errors( + SimulationErrors::PHASEGUESS_TOO_LARGE, temp.netlistInfo.label_);*/ } } // (hbar / 2 * e) ( -(2 / h) φp1 + (1 / 2h) φp2 ) if(atyp_ == AnalysisType::Voltage) { b_.at(temp.variableIndex_) = - (Constants::SIGMA) * (-(2.0 / (stepSize_ * factor)) * - temp.pn1_ + (1.0 / (2.0 * (stepSize_ * factor))) * temp.pn2_); + (Constants::SIGMA) * (-(2.0 / (stepSize_)) * + temp.pn1_ + (1.0 / (2.0 * (stepSize_))) * temp.pn2_); // (4 / 3) φp1 - (1/3) φp2 } else if (atyp_ == AnalysisType::Phase) { b_.at(temp.variableIndex_) = @@ -380,26 +392,43 @@ void Simulation::handle_jj( temp.vn4_ = temp.vn3_; temp.vn3_ = temp.vn2_; // Update junction transition - if(model.value().get_resistanceType() == 1) { + if(model.rtype() == 1) { auto testLU = temp.update_value(v0); if(testLU && !needsLU_) { needsLU_ = true; } } - // -(hR / h + 2RC) * (Ic sin φ0 - 2C / h Vp1 + C/2h Vp2 + It) - b_.at(temp.indexInfo.currentIndex_.value()) = - (temp.matrixInfo.nonZeros_.back()) * ((((Constants::PI * temp.del_) / - (2 * Constants::EV * temp.rncalc_)) * (sin(temp.phi0_) / - sqrt(1 - model.value().get_transparency() * (sin(temp.phi0_ / 2) * - sin(temp.phi0_ / 2)))) * tanh((temp.del_) / - (2 * Constants::BOLTZMANN * model.value().get_temperature()) * - sqrt(1 - model.value().get_transparency() * - (sin(temp.phi0_ / 2) * sin(temp.phi0_ / 2))))) - - (((2 * model.value().get_capacitance()) / - (stepSize_ * factor)) * temp.vn1_) + - ((model.value().get_capacitance() / - (2.0 * (stepSize_ * factor))) * temp.vn2_) + - temp.transitionCurrent_); + // sin (φ0 - φ) + double sin_phi = sin(temp.phi0_ - temp.model_.phiOff()); + if (!temp.model_.tDep()) { + // -(hR / h + 2RC) * (Ic sin (φ0) - 2C / h Vp1 + C/2h Vp2 + It) + b_.at(temp.indexInfo.currentIndex_.value()) = + (temp.matrixInfo.nonZeros_.back()) * (temp.model_.ic() * sin_phi - + (((2 * model.c()) / (stepSize_)) * temp.vn1_) + + ((model.c() / (2.0 * (stepSize_))) * temp.vn2_) + + temp.it_); + } else { + double sin2_half_phi = sin((temp.phi0_ - temp.model_.phiOff()) / 2); + sin2_half_phi = sin2_half_phi * sin2_half_phi; + double sqrt_part = sqrt(1 - model.d() * sin2_half_phi); + b_.at(temp.indexInfo.currentIndex_.value()) = + // -(hR / h + 2RC) *( + (temp.matrixInfo.nonZeros_.back()) * (( + // (π * Δ / 2 * e * Rn) + ((Constants::PI * temp.del_) + / (2 * Constants::EV * temp.model_.rn())) + // * (sin(φ0 - φ) / √(1 - D * sin²((φ0 - φ) / 2)) + * (sin_phi / sqrt_part) + // * tanh(Δ / (2 * kB * T) * √(1 - D * sin²((φ0 - φ) / 2))) + * tanh(temp.del_ / (2 * Constants::BOLTZMANN * model.t()) + * sqrt_part)) + // - 2C / h Vp1 + - (((2 * model.c()) / stepSize_) * temp.vn1_) + // + C/2h Vp2 + + ((model.c() / (2.0 * stepSize_)) * temp.vn2_) + // + It) + + temp.it_); + } temp.vn2_ = temp.vn1_; } } diff --git a/src/Transient.cpp b/src/Transient.cpp index 4d47dff..073a3c2 100644 --- a/src/Transient.cpp +++ b/src/Transient.cpp @@ -9,16 +9,21 @@ using namespace JoSIM; -void Transient::identify_simulation( - const std::vector& controls, Transient& tObj) { +void Transient::identify_simulation( + std::vector& controls, Transient& tObj) { // Flag to store that a transient simulation was found bool transFound = false; // Loop through all the controls - for (const auto& i : controls) { + for (auto& i : controls) { // If the first token of the line contains "TRAN" if (i.front().find("TRAN") != std::string::npos) { // Set the flag as true transFound = true; + // Check for disable startup flag + if (i.back() == "DST") { + tObj.startup(false); + i.pop_back(); + } // If there are less than 2 tokens if (i.size() < 2) { // Complain of invalid transient analysis specification @@ -69,11 +74,10 @@ void Transient::identify_simulation( tObj.tstop(1E-9); tObj.prstart(0); } - // While variable time step does not work yet - // If user provided time step is larger than 0.25ps - if (tObj.tstep() > 0.25E-12) { - tObj.tstep(0.25E-12); - } + //// If user provided time step is larger than 0.25ps junction will fail + //if (tObj.tstep() > 0.25E-12) { + // tObj.tstep(0.25E-12); + //} // Also if PSTEP is smaller than TSTEP if (tObj.tstep() > tObj.prstep()) { // Reduce TSTEP to match PRSTEP diff --git a/src/Verbose.cpp b/src/Verbose.cpp index 1e6a240..c1a0e88 100644 --- a/src/Verbose.cpp +++ b/src/Verbose.cpp @@ -44,6 +44,12 @@ void Verbose::print_circuit_stats(const Input& iObj, const Matrix& mObj) { // Print the total component count std::cout << std::left << std::setw(26) << "Component count:" << mObj.components.devices.size() << "\n"; + // Print the total non-zero count + std::cout << std::left << std::setw(26) << "Non-zero count:" << + mObj.nz.size() << "\n"; + // Print matrix dimensions + std::cout << std::left << std::setw(26) << "MxN size:" << + mObj.rp.size() << "x" << mObj.rp.size() << "\n"; // Print the total jj count if (mObj.components.junctionIndices.size() != 0) { std::cout << std::left << std::setw(26) << "JJ count:" << diff --git a/src/josim.cpp b/src/josim.cpp index a3f8978..86242f5 100644 --- a/src/josim.cpp +++ b/src/josim.cpp @@ -47,16 +47,16 @@ int main(int argc, // Identify the simulation parameters Transient::identify_simulation(iObj.controls, iObj.transSim); // Add noise (if any) - Noise::add_noise_sources(iObj); + //Noise::add_noise_sources(iObj); // Create matrix object Matrix mObj; // Create the matrix in csr format mObj.create_matrix(iObj); // Do verbosity Verbose::handle_verbosity(iObj.argVerb, iObj, mObj); - // Dump expanded Netlist since it is no longer needed - iObj.netlist.expNetlist.clear(); - iObj.netlist.expNetlist.shrink_to_fit(); + //// Dump expanded Netlist since it is no longer needed + //iObj.netlist.expNetlist.clear(); + //iObj.netlist.expNetlist.shrink_to_fit(); // Find the relevant traces to store find_relevant_traces(iObj, mObj); // Create a simulation object