Implement modified Newton

Author: daklauss

src/libcadet/SimulatorImpl.cpp

@@ -225,19 +225,32 @@ namespace cadet
 			cadet::AdJacobianParams{sim->_vecADres, sim->_vecADy, sim->numSensitivityAdDirections()});
 	}
 
-	/*
-	int jacobianUpdateWrapper(N_Vector y, N_Vector yDot, N_Vector res, N_Vector tempv1, N_Vector tempv2, N_Vector tempv3, void* userData)
+	/**
+	* @brief IDAS jacobian wrapper function to call the model's jacobian() method
+	*/
+
+	int jacobianUpdateWrapper(SUNLinearSolver LS, SUNMatrix A)
 	{
-		cadet::Simulator* const sim = static_cast<cadet::Simulator*>(userData);
-		const double t = IDA_mem->ida_tn;
+		cadet::Simulator* const sim = static_cast<cadet::Simulator*>(LS->content);
+
 		const unsigned int secIdx = sim->getCurrentSection(t);
+		double t;
+		double alpha;
+		N_Vector y;
+		N_Vector yDot;
+		N_Vector unused1;
+		N_Vector unused2;
+		N_Vector unused3;
+		void* unused4;
+		IDAGetErrWeights(sim->_idaMemBlock, sim->_linearSolverWeight);
+		IDAGetNonlinearSystemData(sim->_idaMemBlock, &t, &unused1, &unused2, &y, &yDot, &unused3, &alpha, &unused4);
 
 		LOG(Trace) << "==> Jacobian at t = " << t;
 
 		return sim->_model->jacobian(cadet::SimulationTime{ t, secIdx }, cadet::ConstSimulationState{ NVEC_DATA(y), NVEC_DATA(yDot) }, NVEC_DATA(tempv1),
 			cadet::AdJacobianParams{ sim->_vecADres, sim->_vecADy, sim->numSensitivityAdDirections() });
 	}
-	*/
+	
 
 	/**
 	* @brief Change the error weights in the state vector
@@ -367,21 +380,6 @@ namespace cadet
 		return 0;
 	}
 
-	int linearSolverSetup(SUNLinearSolver, SUNMatrix)
-	{
-		return 0;
-	}
-
-	int linearSolverSetATimes(SUNLinearSolver, void*, SUNATimesFn)
-	{
-		return 0;
-	}
-
-	int linearSolverNumIters(SUNLinearSolver)
-	{
-		return 1;
-	}
-
 	N_Vector linearSolverResidual(SUNLinearSolver)
 	{
 		return nullptr;
@@ -427,7 +425,7 @@ namespace cadet
 	Simulator::Simulator() : _model(nullptr), _solRecorder(nullptr), _idaMemBlock(nullptr), _vecStateY(nullptr),
 		_vecStateYdot(nullptr), _vecFwdYs(nullptr), _vecFwdYsDot(nullptr),
 		_relTolS(1.0e-9), _absTol(1, 1.0e-12), _relTol(1.0e-9), _initStepSize(1, 1.0e-6), _maxSteps(10000), _maxStepSize(0.0),
-		_nThreads(0), _sensErrorTestEnabled(true), _maxNewtonIter(4), _maxErrorTestFail(10), _maxConvTestFail(10),
+		_nThreads(0), _modifiedNewton(false), _sensErrorTestEnabled(true), _maxNewtonIter(4), _maxErrorTestFail(10), _maxConvTestFail(10),
 		_maxNewtonIterSens(4), _curSec(0), _skipConsistencyStateY(false), _skipConsistencySensitivity(false),
 		_consistentInitMode(ConsistentInitialization::Full), _consistentInitModeSens(ConsistentInitialization::Full),
 		_vecADres(nullptr), _vecADy(nullptr), _lastIntTime(0.0), _notification(nullptr), _linearSolver(nullptr), _sunctx(nullptr)
@@ -533,6 +531,8 @@ namespace cadet
 		_linearSolver->content = this;
 		_linearSolver->ops->gettype = linearSolverGetType;
 		_linearSolver->ops->solve = linearSolverSolve;
+		if (_modifiedNewton)
+			_linearSolver->ops->setup = jacobianUpdateWrapper;
 
 		// Attach user data structure
 		IDASetUserData(_idaMemBlock, this);

src/libcadet/SimulatorImpl.hpp

@@ -36,7 +36,7 @@ int residualDaeWrapper(double t, N_Vector y, N_Vector yDot, N_Vector res, void*
 int linearSolverSolve(SUNLinearSolver ls, SUNMatrix, N_Vector x, N_Vector b, double tol);
 int linearSolverSetScalingVectors(SUNLinearSolver ls, N_Vector weight, N_Vector);
 
-int jacobianUpdateWrapper(void* IDA_mem, N_Vector y, N_Vector yDot, N_Vector res, N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
+int jacobianUpdateWrapper(SUNLinearSolver LS, SUNMatrix A);
 
 int residualSensWrapper(int ns, double t, N_Vector y, N_Vector yDot, N_Vector res, 
 		N_Vector* yS, N_Vector* ySDot, N_Vector* resS,
@@ -219,7 +219,7 @@ class Simulator : public ISimulator
 
 	friend int ::cadet::linearSolverSolve(SUNLinearSolver ls, SUNMatrix, N_Vector x, N_Vector b, double tol);
 	friend int ::cadet::linearSolverSetScalingVectors(SUNLinearSolver ls, N_Vector weight, N_Vector);
-	friend int ::cadet::jacobianUpdateWrapper(void* IDA_mem, N_Vector y, N_Vector yDot, N_Vector res, N_Vector tempv1, N_Vector tempv2, N_Vector tempv3);
+	friend int ::cadet::jacobianUpdateWrapper(SUNLinearSolver ls, SUNMatrix A);
 
 //	friend int ::cadet::weightWrapper(N_Vector y, N_Vector ewt, void *user_data);