- slightly changes to the implicit radau codes

- added initial testset for radau codes git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@15265 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
OpenModelica · Feb 22, 2013 · 9260697 · 9260697
1 parent a7a2b72
commit 9260697
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 15 deletions.
diff --git a/SimulationRuntime/c/simulation/solver/radau.c b/SimulationRuntime/c/simulation/solver/radau.c
@@ -327,13 +327,14 @@ int initKinsol(KINODE *kinOde)
   double* seq,*slow,*sup;
   KDATAODE *kData = kinOde->kData;
   NLPODE *nlp = kinOde->nlp;
-  double*derx = nlp->derx;
+  double* derx;
   double tmp,tmp1,h;
   double* xeq = NV_DATA_S(kData->x);
   nlp->currentStep = &kinOde->solverInfo->currentStepSize;
   kinOde->nlp->derx = data->localData[0]->realVars + nStates;
   kinOde->nlp->x0 = data->localData[1]->realVars;
   kinOde->nlp->t0 = data->localData[1]->timeValue;
+  derx = nlp->derx;
 
   xlow = xeq + N*nStates;
   xup = xlow + N*nStates;
@@ -395,24 +396,24 @@ int initKinsol(KINODE *kinOde)
 
   kData->mset = 1;
 
-  kData->fnormtol = h;
-  kData->scsteptol = h;
+  kData->fnormtol = data->simulationInfo.tolerance;
+  kData->scsteptol = data->simulationInfo.tolerance;
 
   KINSetFuncNormTol(kData->kmem, kData->fnormtol);
   KINSetScaledStepTol(kData->kmem, kData->scsteptol);
 
   return 0;
 }
 
-int refreshModell(DATA* data, double*x, double time)
+int refreshModell(DATA* data, double* x, double time)
 {
   int i;
   SIMULATION_DATA *sData = (SIMULATION_DATA*)data->localData[0];
 
-  for(i=0;i<data->modelData.nStates;i++)
-    sData->realVars[i] = x[i];
+  memcpy(sData->realVars, x, sizeof(double)*data->modelData.nStates);
   sData->timeValue = time;
   functionODE(data);
+
   return 0;
 }
 
@@ -748,6 +749,7 @@ int lobatto6Res(N_Vector x, N_Vector f, void* user_data)
 
 int kinsolOde(void* ode)
 {
+  int retValue;
   KINODE *kinOde = (KINODE*) ode;
   KDATAODE *kData = kinOde->kData;
   NLPODE *nlp = kinOde->nlp;
@@ -757,23 +759,24 @@ int kinsolOde(void* ode)
   initKinsol(kinOde);
   do{
     /* Call KINDense to specify the linear solver */
-    KINDense(kData->kmem, 3*N*nlp->nStates);
+    KINSpbcg(kData->kmem, 3*N*nlp->nStates);
     KINSetMaxSetupCalls(kData->kmem, kData->mset);
 
     kData->error_code = KINSol( kData->kmem,           /* KINSol memory block */
                                 kData->x,              /* initial guess on input; solution vector */
                                 kData->glstr,          /* global stragegy choice */
                                 kData->sVars,          /* scaling vector, for the variable cc */
                                 kData->sEqns );
-  }while(kData->error_code < 0 && (++kData->glstr) <= KIN_LINESEARCH);
+    retValue = kData->error_code;
+  }while(retValue < 0 && (++kData->glstr) <= KIN_LINESEARCH);
   refreshModell(kinOde->data, NV_DATA_S(kData->x) + (N-1)*nlp->nStates, nlp->t0 + (nlp->dt));
-  return 0;
+  return retValue;
 }
 #else
 
 int kinsolOde(void* ode)
 {
   assert(0);
-  return 0;
+  return -1;
 }
 #endif
diff --git a/SimulationRuntime/c/simulation/solver/radau.h b/SimulationRuntime/c/simulation/solver/radau.h
@@ -35,9 +35,9 @@
 #ifndef _RADAU_H_
 #define _RADAU_H_
 
-#include "../../../Compiler/runtime/config.h"
 #include "simulation_data.h"
 #include "solver_main.h"
+#include "../../Compiler/runtime/config.h"
 
 #ifdef WITH_SUNDIALS
   #include <math.h>

diff --git a/SimulationRuntime/c/simulation/solver/solver_main.c b/SimulationRuntime/c/simulation/solver/solver_main.c
@@ -30,7 +30,8 @@
 
  /*! \file solver_main.c
  */
-
+
+#include "../../Compiler/runtime/config.h"
 #include "solver_main.h"
 #include "simulation_runtime.h"
 #include "simulation_result.h"
@@ -726,12 +727,15 @@ int solver_main(DATA* data, const char* init_initMethod,
   if (!retVal)
     retVal = initializeModel(data, init_initMethod, init_optiMethod, init_file, init_time, lambda_steps);
 
-  INFO(LOG_SOLVER, "Performed initial value calculation.");
-  INFO2(LOG_SOLVER, "Start numerical solver from %g to %g", simInfo->startTime, simInfo->stopTime);
-
   /* starts the simulation main loop */
   if (!retVal)
+  {
+    INFO(LOG_SOLVER, "Performed initialization.");
+    INFO2(LOG_SOLVER, "Start numerical solver from %g to %g", simInfo->startTime, simInfo->stopTime);
     retVal = performSimulation(data, &solverInfo);
+  }
+  else
+    WARNING(LOG_STDOUT, "Initialization failed.");
 
   /* terminate the simulation */
   finishSimulation(data, &solverInfo, outputVariablesAtEnd);