Skip to content

Commit

Permalink
- radau, lobatto code style
Browse files Browse the repository at this point in the history 
git-svn-id: https://openmodelica.org/svn/OpenModelica/trunk@15414 f25d12d1-65f4-0310-ae8a-bbce733d8d8e
  • Loading branch information
Vitalij Ruge committed Mar 1, 2013
1 parent 117822f commit 994fb6a
Show file tree
Hide file tree
Showing 4 changed files with 48 additions and 50 deletions.
4 changes: 4 additions & 0 deletions SimulationRuntime/c/simulation/simulation_runtime.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -735,7 +735,11 @@ int callSolver(DATA* simData, string result_file_cstr, string init_initMethod,
#endif
} else {
INFO1(LOG_STDOUT, " | Unrecognized solver: %s.", simData->simulationInfo.solverMethod);
#ifdef WITH_SUNDIALS
INFO(LOG_STDOUT, " | valid solvers are: dassl, euler, rungekutta, inline-euler, inline-rungekutta, dasslwort, dasslSymJac, dasslNumJac, dasslColorSymJac, dasslInternalNumJac, qss, radau1, radau3, radau5, lobatto2, lobatto4");
#else
INFO(LOG_STDOUT, " | valid solvers are: dassl, euler, rungekutta, inline-euler, inline-rungekutta, dasslwort, dasslSymJac, dasslNumJac, dasslColorSymJac, dasslInternalNumJac, qss.");
#endif
retVal = 1;
}

Expand Down
79 changes: 37 additions & 42 deletions SimulationRuntime/c/simulation/solver/radau.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -251,6 +251,7 @@ static int allocateKINSOLODE(KINODE *kinOde)
int nn = 3*N*n; /* 3*N = N*(subintervalls) * 3*var(eq,low,up) */
int i, j, k;
double* ceq,*clow,*cup;
double* seq,*slow,*sup;
double* sveq,*svlow,*svup;
KDATAODE * kData = (kinOde)->kData;
NLPODE *nlp = kinOde->nlp;
Expand All @@ -266,6 +267,9 @@ static int allocateKINSOLODE(KINODE *kinOde)
sveq = NV_DATA_S(kData->sVars);
svlow = sveq + N*n;
svup = svlow + N*n;
seq = NV_DATA_S(kData->sEqns);
slow = seq + N*n;
sup = slow + N*n;

for(j=0, k=0; j<N; j++)
{
Expand All @@ -277,6 +281,8 @@ static int allocateKINSOLODE(KINODE *kinOde)
sveq[k] = s[i];
svlow[k] = s[i];
svup[k] = s[i];
slow[k] = s[i];
sup[k] = s[i];
}
}
KINSetUserData(kinOde->kData->kmem, (void*) kinOde);
Expand Down Expand Up @@ -395,8 +401,6 @@ static int initKinsol(KINODE *kinOde)
xup[k] = xeq[k] - nlp->max[i];
tmp = 1.0/(fabs(nlp->x0[i] - xeq[k]) + 1e-6);
seq[k] = tmp;
slow[k] = tmp;
sup[k] = tmp;
}
}

Expand Down Expand Up @@ -425,7 +429,7 @@ static int refreshModell(DATA* data, double* x, double time)

static int radau5Res(N_Vector x, N_Vector f, void* user_data)
{
int i,sub2,sub3;
int i,k;
KINODE* kinOde = (KINODE*)user_data;
NLPODE *nlp = kinOde->nlp;
DATA *data = kinOde->data;
Expand Down Expand Up @@ -463,40 +467,36 @@ static int radau5Res(N_Vector x, N_Vector f, void* user_data)
feq[i] = (nlp->c[0][0]*x0[i] + nlp->c[0][3]*x3[i] + h*derx[i]) -
(nlp->c[0][1]*x1[i] + nlp->c[0][2]*x2[i]);


flow[i] = xlow[i] - x1[i] + lb[i];
fup[i] = xup[i] - x1[i] + ub[i];
}

refreshModell(data, x2,t0 + a[1]*h);
for(i = 0;i<n;i++)
for(i = 0, k=n; i<n; i++, k++)
{
sub2 = i + n;

feq[sub2] = (nlp->c[1][1]*x1[i] + h*derx[i]) -
feq[k] = (nlp->c[1][1]*x1[i] + h*derx[i]) -
(nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i] + nlp->c[1][3]*x3[i]);

flow[sub2] = xlow[sub2] - x2[i] + lb[i];
fup[sub2] = xup[sub2] - x2[i] + ub[i];
flow[k] = xlow[k] - x2[i] + lb[i];
fup[k] = xup[k] - x2[i] + ub[i];
}

refreshModell(data, x3, t0 + a[2]*h);
for(i = 0;i<n;i++)
refreshModell(data, x3, t0 + h);
for(i = 0;i<n;i++,k++)
{
sub3 = i + 2*n;
feq[sub3] = (nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i] + h*derx[i]) -
feq[k] = (nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i] + h*derx[i]) -
(nlp->c[2][1]*x1[i] + nlp->c[2][3]*x3[i]);

flow[sub3] = xlow[sub3] - x3[i] + lb[i];
fup[sub3] = xup[sub3] - x3[i] + ub[i];
flow[k] = xlow[k] - x3[i] + lb[i];
fup[k] = xup[k] - x3[i] + ub[i];
}

return 0;
}

static int radau3Res(N_Vector x, N_Vector f, void* user_data)
{
int i,sub2, N;
int i,k, N;
KINODE* kinOde = (KINODE*)user_data;
NLPODE *nlp = kinOde->nlp;
DATA *data = kinOde->data;
Expand Down Expand Up @@ -534,15 +534,13 @@ static int radau3Res(N_Vector x, N_Vector f, void* user_data)
fup[i] = xup[i] - x1[i] + ub[i];
}

refreshModell(data, x2,t0 + a[1]*h);
for(i = 0;i<n;i++)
refreshModell(data, x2,t0 + h);
for(i = 0, k=n;i<n;i++,k++)
{
sub2 = i + n;

feq[sub2] = (nlp->c[1][1]*x1[i] + h*derx[i]) -
feq[k] = (nlp->c[1][1]*x1[i] + h*derx[i]) -
(nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i]);
flow[sub2] = xlow[sub2] - x2[i] + lb[i];
fup[sub2] = xup[sub2] - x2[i] + ub[i];
flow[k] = xlow[k] - x2[i] + lb[i];
fup[k] = xup[k] - x2[i] + ub[i];
}

return 0;
Expand All @@ -551,7 +549,7 @@ static int radau3Res(N_Vector x, N_Vector f, void* user_data)

static int radau1Res(N_Vector x, N_Vector f, void* user_data)
{
int i,sub2, N;
int i, N;
KINODE* kinOde = (KINODE*)user_data;
NLPODE *nlp = kinOde->nlp;
DATA *data = kinOde->data;
Expand Down Expand Up @@ -631,7 +629,7 @@ static int lobatto2Res(N_Vector x, N_Vector f, void* user_data)

static int lobatto4Res(N_Vector x, N_Vector f, void* user_data)
{
int i,sub2, N;
int i,k, N;
KINODE* kinOde = (KINODE*)user_data;
NLPODE *nlp = kinOde->nlp;
DATA *data = kinOde->data;
Expand Down Expand Up @@ -671,19 +669,18 @@ static int lobatto4Res(N_Vector x, N_Vector f, void* user_data)
}

refreshModell(data, x2,t0 + h);
for(i = 0;i<n;i++)
for(i = 0,k=n;i<n;i++,k++)
{
sub2 = i + n;
feq[sub2] = (2.0*h*derx[i] + 16.0*x1[i]) - (8.0*(x0[i] + x2[i]) +2.0*h*f0[i]);
flow[sub2] = xlow[sub2] - x2[i] + lb[i];
fup[sub2] = xup[sub2] - x2[i] + ub[i];
feq[k] = (2.0*h*derx[i] + 16.0*x1[i]) - (8.0*(x0[i] + x2[i]) +2.0*h*f0[i]);
flow[k] = xlow[k] - x2[i] + lb[i];
fup[k] = xup[k] - x2[i] + ub[i];
}
return 0;
}

static int lobatto6Res(N_Vector x, N_Vector f, void* user_data)
{
int i,sub2, N;
int i,k, N;
KINODE* kinOde = (KINODE*)user_data;
NLPODE *nlp = kinOde->nlp;
DATA *data = kinOde->data;
Expand Down Expand Up @@ -725,21 +722,19 @@ static int lobatto6Res(N_Vector x, N_Vector f, void* user_data)
}

refreshModell(data, x2,t0 + nlp->a[1]*h);
for(i = 0;i<n;i++)
for(i = 0,k=n;i<n;i++,k++)
{
sub2 = i + n;
feq[sub2] = (h*derx[i] + nlp->c[1][1]*x1[i]) - (h*nlp->c[1][4]*f0[i] + nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i] + nlp->c[1][3]*x3[i]);
flow[sub2] = xlow[sub2] - x2[i] + lb[i];
fup[sub2] = xup[sub2] - x2[i] + ub[i];
feq[k] = (h*derx[i] + nlp->c[1][1]*x1[i]) - (h*nlp->c[1][4]*f0[i] + nlp->c[1][0]*x0[i] + nlp->c[1][2]*x2[i] + nlp->c[1][3]*x3[i]);
flow[k] = xlow[k] - x2[i] + lb[i];
fup[k] = xup[k] - x2[i] + ub[i];
}

refreshModell(data, x3,t0 + h);
for(i = 0;i<n;i++)
for(i = 0;i<n;i++,k++)
{
sub2 = i + 2*n;
feq[sub2] = (h*(f0[i] + derx[i]) + nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i]) - (nlp->c[2][1]*x1[i] + nlp->c[2][3]*x3[i]);
flow[sub2] = xlow[sub2] - x3[i] + lb[i];
fup[sub2] = xup[sub2] - x3[i] + ub[i];
feq[k] = (h*(f0[i] + derx[i]) + nlp->c[2][0]*x0[i] + nlp->c[2][2]*x2[i]) - (nlp->c[2][1]*x1[i] + nlp->c[2][3]*x3[i]);
flow[k] = xlow[k] - x3[i] + lb[i];
fup[k] = xup[k] - x3[i] + ub[i];
}

return 0;
Expand Down
2 changes: 1 addition & 1 deletion SimulationRuntime/c/simulation/solver/radau.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,7 +104,7 @@
#endif /* SUNDIALS */
int allocateKinOde(DATA* data, SOLVER_INFO* solverInfo, int flag, int N);
int freeKinOde(DATA* data, SOLVER_INFO* solverInfo, int flag, int N);

int kinsolOde(void* ode);
#ifdef __cplusplus
};
#endif
Expand Down
13 changes: 6 additions & 7 deletions SimulationRuntime/c/simulation/solver/solver_main.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -815,48 +815,47 @@ static int rungekutta_step(DATA* data, SOLVER_INFO* solverInfo)
/*************************************** Radau5 IIA ***********************************/
int radau5_step(DATA* data, SOLVER_INFO* solverInfo)
{

kinsolOde(solverInfo->solverData, 3, data, 6);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}

/*************************************** Radau3 IIA ***********************************/
int radau3_step(DATA* data, SOLVER_INFO* solverInfo)
{
kinsolOde(solverInfo->solverData, 2, data, 7);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}

/*************************************** Radau1 IIA ***********************************/
int radau1_step(DATA* data, SOLVER_INFO* solverInfo)
{
kinsolOde(solverInfo->solverData, 1, data, 8);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}

/*************************************** Lobatto2 IIA ***********************************/
int lobatto2_step(DATA* data, SOLVER_INFO* solverInfo)
{
kinsolOde(solverInfo->solverData, 1, data, 9);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}

/*************************************** Lobatto4 IIA ***********************************/
int lobatto4_step(DATA* data, SOLVER_INFO* solverInfo)
{
kinsolOde(solverInfo->solverData, 2, data, 10);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}

/*************************************** Lobatto6 IIA ***********************************/
int lobatto6_step(DATA* data, SOLVER_INFO* solverInfo)
{
kinsolOde(solverInfo->solverData, 3, data, 11);
kinsolOde(solverInfo->solverData);
solverInfo->currentTime += solverInfo->currentStepSize;
return 0;
}
Expand Down

0 comments on commit 994fb6a

Please sign in to comment.