diff --git a/src/specificClass.js b/src/specificClass.js
index fc31fdc..5f6e2c7 100644
--- a/src/specificClass.js
+++ b/src/specificClass.js
@@ -10,9 +10,9 @@ const mathConfig = {
 
 const math = create(all, mathConfig);
 
-const BC_PADDING = 2;
+const BC_PADDING = 2; // Boundary Condition padding for open boundaries in extendedState variable
 const DEBUG = false;
-const DAY2MS = 1000 * 60 * 60 * 24;
+const DAY2MS = 1000 * 60 * 60 * 24; // convert between days and milliseconds
 
 class Reactor {
   /**
@@ -25,13 +25,14 @@ class Reactor {
     this.logger = new logger(this.config.general.logging.enabled, this.config.general.logging.logLevel, config.general.name);
     this.emitter = new EventEmitter();
     this.measurements = new MeasurementContainer();
-    this.childRegistrationUtils = new childRegistrationUtils(this); // Child registration utility
+    this.childRegistrationUtils = new childRegistrationUtils(this); // child registration utility
 
+    // placeholder variables for children and parents
     this.upstreamReactor = null;
     this.downstreamReactor = null;
     this.returnPump = null;
 
-    this.asm = new ASM3();
+    this.asm = new ASM3(); // Reaction model
 
     this.volume = config.volume; // fluid volume reactor [m3]
 
@@ -44,7 +45,7 @@ class Reactor {
 
     this.currentTime = null; // milliseconds since epoch [ms]
     this.timeStep = 1 / (24*60*60) * this.config.timeStep; // time step in seconds, converted to days.
-    this.speedUpFactor = 100; // speed up factor for simulation, 60 means 1 minute per simulated second
+    this.speedUpFactor = 1; // speed up factor for simulation, 60 means 1 minute per simulated second
   }
 
   /**
@@ -113,6 +114,11 @@ class Reactor {
     }
   }
 
+  /**
+   * Register child function required for child registration.
+   * @param {object} child 
+   * @param {string} softwareType 
+   */
   registerChild(child, softwareType) {
     if(!child) {
       this.logger.error(`Invalid ${softwareType} child provided.`);
@@ -161,14 +167,14 @@ class Reactor {
 
   _connectReactor(reactorChild) {
     if (reactorChild.config.functionality.positionVsParent != "upstream") {
-      this.logger.warn("Reactor children of reactors should always be upstream.");
+      this.logger.warn("Reactor children of other reactors should always be upstream!");
     }
 
     // set upstream and downstream reactor variable in current and child nodes respectively for easy access
     this.upstreamReactor = reactorChild;
     reactorChild.downstreamReactor = this;
 
-    reactorChild.emitter.on("stateChange", (eventData) => {
+    reactorChild.emitter.on("stateChange", (eventData) => { // Triggers state update in downstream reactor.
       this.logger.debug(`State change of upstream reactor detected.`);
       this.updateState(eventData);
     });
@@ -199,21 +205,19 @@ class Reactor {
    * Update the reactor state based on the new time.
    * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
    */
-  updateState(newTime) { // expect update with timestamp
-    if (!this.currentTime) {
+  updateState(newTime) {
+    if (!this.currentTime) { // initialise currentTime variable
       this.currentTime = newTime;
       return;
     }
 
-    if (this.upstreamReactor) {
-      // grab main effluent upstream reactor
+    if (this.upstreamReactor) { // grab main effluent upstream reactor
       this.setInfluent = this.upstreamReactor.getEffluent[0];
     }
 
     const n_iter = Math.floor(this.speedUpFactor * (newTime-this.currentTime) / (this.timeStep*DAY2MS));
 
-    if (n_iter == 0) {
-      // no update required
+    if (n_iter == 0) { // no update required, change in currentTime smaller than time step
       return;
     }
 
@@ -223,9 +227,9 @@ class Reactor {
       n += 1;
     }
     this.currentTime += n_iter * this.timeStep * DAY2MS / this.speedUpFactor;
-    this.emitter.emit("stateChange", this.currentTime);
+    this.emitter.emit("stateChange", this.currentTime); // update downstream reactors
 
-    if (this.returnPump) {
+    if (this.returnPump) { // update recirculation pump state
       this.returnPump.updateSourceSink();
     }
   }
@@ -324,13 +328,16 @@ class Reactor_PFR extends Reactor {
     super._connectReactor(reactorChild);
   }
 
+  /**
+   * Update the reactor state based on the new time. Performs checks specific to PFR.
+   * @param {number} newTime - New time to update reactor state to, in milliseconds since epoch.
+   */
   updateState(newTime) {
     super.updateState(newTime);
-    // let Pe_local = this.d_x*math.sum(this.Fs)/(this.D*this.A)
-    this.D = this._constrainDispersion(this.D);
+    
+    this.D = this._constrainDispersion(this.D); // constrains D to minimum dispersion, so that local Péclet number is always above 2
     const Co_D = this.D*this.timeStep/(this.d_x*this.d_x);
 
-    // (Pe_local >= 2) && this.logger.warn(`Local Péclet number (${Pe_local}) is too high! Increase reactor resolution.`);
     (Co_D >= 0.5) && this.logger.warn(`Courant number (${Co_D}) is too high! Reduce time step size.`);
 
     if(DEBUG) {
@@ -409,8 +416,8 @@ class Reactor_PFR extends Reactor {
    */
   _applyBoundaryConditions() {
     // Upstream BC
-    if (this.upstreamReactor) {
-      // Open boundary
+    if (this.upstreamReactor && this.upstreamReactor.config.reactor_type == "PFR") {
+      // Open boundary, if upstream reactor is PFR
       this.extendedState.splice(0, BC_PADDING, ...this.upstreamReactor.state.slice(-BC_PADDING)); 
     } else {
       if (math.sum(this.Fs) > 0) {
@@ -418,7 +425,7 @@ class Reactor_PFR extends Reactor {
         const BC_C_in = math.multiply(1 / math.sum(this.Fs), [this.Fs], this.Cs_in)[0];
         const BC_dispersion_term = this.D*this.A/(math.sum(this.Fs)*this.d_x);
         this.extendedState[BC_PADDING] = math.multiply(1/(1+BC_dispersion_term), math.add(BC_C_in, math.multiply(BC_dispersion_term, this.extendedState[BC_PADDING+1])));
-        // Numerical boundary condition
+        // Numerical boundary condition (first-order accurate)
         this.extendedState[BC_PADDING-1] = math.add(math.multiply(2, this.extendedState[BC_PADDING]), math.multiply(-2, this.extendedState[BC_PADDING+2]), this.extendedState[BC_PADDING+3]);
       } else {
         // Neumann BC (no flux)
@@ -427,8 +434,8 @@ class Reactor_PFR extends Reactor {
     }
 
     // Downstream BC
-    if (this.downstreamReactor) {
-      // Open boundary
+    if (this.downstreamReactor && this.downstreamReactor.config.reactor_type == "PFR") {
+      // Open boundary, if downstream reactor is PFR
       this.extendedState.splice(this.n_x+BC_PADDING, BC_PADDING, ...this.downstreamReactor.state.slice(0, BC_PADDING));
     } else {
       // Neumann BC (no flux)
@@ -438,7 +445,6 @@ class Reactor_PFR extends Reactor {
 
   /**
    * Create finite difference first derivative operator.
-   * @returns {Array} - First derivative operator matrix.
    */
   _makeDoperator() { // create gradient operator
     const D_size = this.n_x+2*BC_PADDING;
@@ -454,7 +460,6 @@ class Reactor_PFR extends Reactor {
 
   /**
    * Create central finite difference second derivative operator.
-   * @returns {Array} - Second derivative operator matrix.
    */
   _makeD2operator() { // create the central second derivative operator
     const D_size = this.n_x+2*BC_PADDING;
@@ -467,6 +472,9 @@ class Reactor_PFR extends Reactor {
     return D2;
   }
 
+  /**
+   * Constrains dispersion so that local Péclet number stays below 2. Needed for stable central differencing method.
+   */
   _constrainDispersion(D) {
     const Dmin = math.sum(this.Fs) * this.d_x / (1.999 * this.A);
     if (D < Dmin) {