CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM

Using ODE in OpenFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   October 9, 2007, 18:00
Default Hi, I am new to OpenFoam an
  #1
Senior Member
 
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 8
skabilan is on a distinguished road
Hi,

I am new to OpenFoam and was trying to use the ODE
solver. I took a look at the ODETest but was not able to understand how the equations are cast. can some one show me how I can solve a simple system like the following?

C --- INITIAL VALUES
X=0.0D0
Y(1)=2.0D0
Y(2)=-0.66D0
C --- ENDPOINT OF INTEGRATION
XEND=2.0D0
XDELTA=0.2
C --- SYSTEM OF EQUATIONS
F(1)=Y(2)
F(2)=((1-Y(1)**2)*Y(2)-Y(1))/2.0


ODETest
\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "IOmanip.H"
#include "ODE.H"
#include "ODESolver.H"
#include "RK.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

class testODE
:
public ODE
{

public:

testODE()
{}

label nEqns() const
{
return 4;
}

void derivatives
(
const scalar x,
const scalarField& y,
scalarField& dydx
) const
{
dydx[0] = -y[1];
dydx[1] = y[0] - (1.0/x)*y[1];
dydx[2] = y[1] - (2.0/x)*y[2];
dydx[3] = y[2] - (3.0/x)*y[3];
}

void jacobian
(
const scalar x,
const scalarField& y,
scalarField& dfdx,
Matrix<scalar>& dfdy
) const
{
dfdx[0] = 0.0;
dfdx[1] = (1.0/sqr(x))*y[1];
dfdx[2] = (2.0/sqr(x))*y[2];
dfdx[3] = (3.0/sqr(x))*y[3];

dfdy[0][0] = 0.0;
dfdy[0][1] = -1.0;
dfdy[0][2] = 0.0;
dfdy[0][3] = 0.0;

dfdy[1][0] = 1.0;
dfdy[1][1] = -1.0/x;
dfdy[1][2] = 0.0;
dfdy[1][3] = 0.0;

dfdy[2][0] = 0.0;
dfdy[2][1] = 1.0;
dfdy[2][2] = -2.0/x;
dfdy[2][3] = 0.0;

dfdy[3][0] = 0.0;
dfdy[3][1] = 0.0;
dfdy[3][2] = 1.0;
dfdy[3][3] = -3.0/x;
}
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
argList::validArgs.clear();
argList::validArgs.append("ODESolver");
argList args(argc, argv);

word ODESolverName(args.args()[1]);

testODE ode;
autoPtr<odesolver> odeSolver = ODESolver::New(ODESolverName, ode);

scalar xStart = 1.0;
scalarField yStart(ode.nEqns());
yStart[0] = ::Foam::j0(xStart);
yStart[1] = ::Foam::j1(xStart);
yStart[2] = ::Foam::jn(2, xStart);
yStart[3] = ::Foam::jn(3, xStart);

scalarField dyStart(ode.nEqns());
ode.derivatives(xStart, yStart, dyStart);

Info<< setw(10) << "eps" << setw(12) << "hEst";
Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl;
Info<< setprecision(6);

for (label i=0; i<15; i++)
{
scalar eps = ::Foam::exp(-scalar(i + 1));

scalar x = xStart;
scalarField y = yStart;
scalarField dydx = dyStart;

scalarField yScale(ode.nEqns(), 1.0);
scalar hEst = 0.6;
scalar hDid, hNext;
odeSolver->solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext);

Info<< scientific << setw(13) << eps;
Info<< fixed << setw(11) << hEst;
Info<< setw(13) << hDid << setw(13) << hNext
<< setw(13) << y[0] << setw(13) << y[1]
<< setw(13) << y[2] << setw(13) << y[3]
<< endl;
}

scalar x = xStart;
scalar xEnd = x + 1.0;
scalarField y = yStart;

scalarField yEnd(ode.nEqns());
yEnd[0] = ::Foam::j0(xEnd);
yEnd[1] = ::Foam::j1(xEnd);
yEnd[2] = ::Foam::jn(2, xEnd);
yEnd[3] = ::Foam::jn(3, xEnd);

scalar hEst = 0.5;

odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst);
Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
Info << "Numerical: y(2.0) = " << y << ", hEst = " << hEst << endl;

Info << "\nEnd\n" << endl;

return 0;
}


// ************************************************** *********************** //
skabilan is offline   Reply With Quote

Old   October 19, 2007, 14:23
Default Hi, I figured out how to so
  #2
Senior Member
 
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 8
skabilan is on a distinguished road
Hi,

I figured out how to solve ODE's using the above driver.

1) Change the dimension of the system.

abel nEqns() const
{
return 4;
}


2) Change the following equations accordingly

{
dydx[0] = -y[1];
dydx[1] = y[0] - (1.0/x)*y[1];
dydx[2] = y[1] - (2.0/x)*y[2];
dydx[3] = y[2] - (3.0/x)*y[3];
}

3) Specifiy the " scalar xStart " , " scalar xEnd " and you can specify the initial values for yStart's if required.

Feel free to contact if not clear!
skabilan is offline   Reply With Quote

Old   October 19, 2007, 14:56
Default Good work! Does the solution c
  #3
Senior Member
 
Srinath Madhavan (a.k.a pUl|)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 698
Rep Power: 12
msrinath80 is on a distinguished road
Good work! Does the solution compare well to the analytical result?
msrinath80 is offline   Reply With Quote

Old   February 19, 2008, 18:01
Default Hello there Does anyone kno
  #4
New Member
 
Bernardo D Flores
Join Date: Mar 2009
Posts: 21
Rep Power: 8
david_flo1 is on a distinguished road
Hello there

Does anyone knows how to integrate the ODE on OpenFOAM, i.e., using a solver like icoFoam and a ODE system serving as input on icoFOAM inlets to be solved at each time step.

Any pointers are greatly appreciated

Regards
david_flo1 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
OpenFoam vs CFX5 mass balance in OpenFoam tangd OpenFOAM Running, Solving & CFD 33 May 23, 2010 16:36
CheckMesh error using a tutorial from OpenFOAM 114 with openFOAM 13 martapajon OpenFOAM Native Meshers: blockMesh 7 January 21, 2008 13:52
OpenFOAM users in Munich OpenFOAM benutzer in M%c3%bcnchen jaswi OpenFOAM 0 August 3, 2007 13:11
A new Howto on the OpenFOAM Wiki Compiling OpenFOAM under Unix mbeaudoin OpenFOAM Installation 2 April 28, 2006 08:54
Is GUI there for OpenFoam navaladi OpenFOAM 3 March 28, 2006 09:08


All times are GMT -4. The time now is 19:30.