Add write last sensitivity state for units

cadet · Apr 19, 2024 · 40682b4 · 40682b4
1 parent 94376e7
commit 40682b4
Show file tree

Hide file tree

Showing 2 changed files with 83 additions and 1 deletion.
diff --git a/doc/interface/return_data.rst b/doc/interface/return_data.rst
@@ -297,4 +297,12 @@ Group /input/return/unit_XXX
    =============  ==========================
    **Type:** int  **Range:** :math:`\{0,1\}`
    =============  ==========================
+
+``WRITE_SENS_LAST_UNIT``
+
+   Write sensitivity state vector of this unit at last time point, including the state derivative vector (optional, defaults to 0)
+
+   =============  ==========================
+   **Type:** int  **Range:** :math:`\{0,1\}`
+   =============  ==========================
 
diff --git a/include/common/Driver.hpp b/include/common/Driver.hpp
@@ -455,7 +455,6 @@ class Driver
 		else
 			_writeLastState = false;
 
-		// Adding a new output variable
 		std::ostringstream oss;
 		for (int i = 0; i <= _sim->model()->maxUnitOperationId(); ++i)
 		{
@@ -481,6 +480,25 @@ class Driver
 		else
 			_writeLastStateSens = false;
 
+		for (int i = 0; i <= _sim->model()->maxUnitOperationId(); ++i)
+		{
+			oss.str("");
+			oss << "unit_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << i;
+			if (pp.exists(oss.str()))
+			{
+				pp.pushScope(oss.str());
+
+				if (pp.exists("WRITE_SENS_LAST_UNIT"))
+				{
+					const bool writeLastStateSensUnit = pp.getBool("WRITE_SENS_LAST_UNIT");
+					if (writeLastStateSensUnit)
+						_writeLastStateSensUnitId.push_back(i);
+				}
+
+				pp.popScope();
+			}
+		}
+
 		pp.popScope(); // scope return
 
 		if (applyInSimulator)
@@ -588,6 +606,40 @@ class Driver
 			}
 		}
 
+		for (int i = 0; i < _writeLastStateSensUnitId.size(); ++i)
+		{
+			unsigned int sliceStart;
+			unsigned int sliceEnd;
+			std::tie(sliceStart, sliceEnd) = _sim->model()->getModelStateOffsets(_writeLastStateSensUnitId[i]);
+
+			oss.str("");
+			oss << "unit_" << std::setfill('0') << std::setw(3) << std::setprecision(0) << _writeLastStateSensUnitId[i];
+			writer.pushGroup("sensitivity");
+			std::string unitID = oss.str();
+
+			unsigned int len = 0;
+
+			const std::vector<double const*> lastSens = _sim->getLastSensitivities(len);
+			const std::vector<double const*> lastSensdot = _sim->getLastSensitivityDerivatives(len);
+
+			const int nSens = lastSens.size();
+
+			for (int i = 0; i < nSens; i++)
+			{
+				oss.str("");
+				oss << std::setfill('0') << std::setw(3) << std::setprecision(0) << i;
+				writer.pushGroup("param_" + oss.str());
+				writer.pushGroup(unitID);
+				writer.template vector<double>("LAST_SENS_Y", sliceEnd - sliceStart, lastSens[i] + sliceStart);
+				writer.template vector<double>("LAST_SENS_YDOT", sliceEnd - sliceStart, lastSensdot[i] + sliceStart);
+				writer.popGroup();
+				writer.popGroup();
+			}
+
+			writer.popGroup();
+			writer.popGroup();
+		}
+
 		writer.popGroup();
 
 		if (writer.exists("meta"))
@@ -651,6 +703,27 @@ class Driver
 		}
 	}
 
+	inline void setWriteLastSensStateOfUnit(UnitOpIdx uid, bool writeLastSensStateUnit)
+	{
+		const std::vector<UnitOpIdx>::iterator it = std::find(_writeLastStateSensUnitId.begin(), _writeLastStateSensUnitId.end(), uid);
+		if (writeLastSensStateUnit)
+		{
+			if (it == _writeLastStateSensUnitId.end())
+			{
+				// Unit not in list, so add it
+				_writeLastStateSensUnitId.push_back(uid);
+			}
+		}
+		else
+		{
+			if (it != _writeLastStateSensUnitId.end())
+			{
+				// Unit is listed, so remove it
+				_writeLastStateSensUnitId.erase(it);
+			}
+		}
+	}
+
 	inline void setWriteLastState(bool writeLastState) CADET_NOEXCEPT { _writeLastState = writeLastState; }
 	inline void setWriteLastStateSens(bool writeLastState) CADET_NOEXCEPT { _writeLastStateSens = writeLastState; }
 	inline void setWriteSolutionTimes(bool solTimes) CADET_NOEXCEPT
@@ -670,6 +743,7 @@ class Driver
 	bool _writeLastState;
 	std::vector<UnitOpIdx> _writeLastStateUnitId;
 	bool _writeLastStateSens;
+	std::vector<UnitOpIdx> _writeLastStateSensUnitId;
 
 	/**
 	 * @brief Sets section times and section continuity from the given parameter provider