# Using ODE in OpenFoam

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

 October 9, 2007, 18:00 Hi, I am new to OpenFoam an #1 Senior Member   Senthil Kabilan Join Date: Mar 2009 Posts: 113 Rep Power: 8 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& 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::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; } // ************************************************** *********************** //

 October 19, 2007, 14:23 Hi, I figured out how to so #2 Senior Member   Senthil Kabilan Join Date: Mar 2009 Posts: 113 Rep Power: 8 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!

 October 19, 2007, 14:56 Good work! Does the solution c #3 Senior Member   Srinath Madhavan (a.k.a pUl|) Join Date: Mar 2009 Location: Edmonton, AB, Canada Posts: 700 Rep Power: 12 Good work! Does the solution compare well to the analytical result?

 February 19, 2008, 18:01 Hello there Does anyone kno #4 New Member   Bernardo D Flores Join Date: Mar 2009 Posts: 21 Rep Power: 8 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

 Yesterday, 19:30 #5 New Member   methma Rajamuni Join Date: Jul 2015 Posts: 4 Rep Power: 2 Hi, Can you please explain in more details the steps of solving an ODE in OpenFOAM? Thanks, Methma

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post tangd OpenFOAM Running, Solving & CFD 33 May 23, 2010 16:36 martapajon OpenFOAM Native Meshers: blockMesh 7 January 21, 2008 13:52 jaswi OpenFOAM 0 August 3, 2007 13:11 mbeaudoin OpenFOAM Installation 2 April 28, 2006 08:54 navaladi OpenFOAM 3 March 28, 2006 09:08

All times are GMT -4. The time now is 22:58.