CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Using ODE in OpenFoam

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

Like Tree10Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   October 9, 2007, 19:00
Default Using ODE in OpenFoam
  #1
Senior Member
 
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 17
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;
}


// ************************************************** *********************** //
ahmmedshakil and Zhiheng Wang like this.
skabilan is offline   Reply With Quote

Old   October 19, 2007, 15:23
Default Hi, I figured out how to so
  #2
Senior Member
 
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 17
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!
Zhiheng Wang likes this.
skabilan is offline   Reply With Quote

Old   October 19, 2007, 15: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: 703
Rep Power: 21
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: 17
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

Old   September 2, 2015, 20:30
Default
  #5
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hi,

Can you please explain in more details the steps of solving an ODE in OpenFOAM?

Thanks,

Methma
meth is offline   Reply With Quote

Old   September 9, 2015, 08:30
Default
  #6
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
Did you try ODE test?

Also check this tutorial
meth likes this.
hk318i is offline   Reply With Quote

Old   September 14, 2015, 03:54
Default
  #7
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hi Hassan ,

I went through the tutorial you mentioned. I changed the ODETest.C file so that it contains my equations. After that I run the wmake command to complete the C file. It was compiled by making dependency list for source file. But I don't know how to get the results of the ODE solver. Can you please help me? I am new to OpenFOAM. I really need to understand this.

My final goal is to work on fluid structure interface (Study on vortex induced vibration of a sphere/cylinder on free stream). There I need to solve the structure motion ODE in each time step.

Thank you.

Methma
meth is offline   Reply With Quote

Old   September 14, 2015, 04:12
Default
  #8
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hi Hassan ,

I went through the tutorial you mentioned. I changed the TestODE.C file so that it contains my equations. After that I run the wmake command to complete the C file. It was compiled by making dependency list for source file. But I don't know how to get the results of the ODE solver. When I compile the C file does it run automatically or do I need to run it separately? Can you please help me? I am new to OpenFOAM. I really need to understand this.

My final goal is to work on fluid structure interface (Study on vortex induced vibration of a sphere/cylinder on free stream). There I need to solve the structure motion ODE in each time step.

Thank you.

Methma
meth is offline   Reply With Quote

Old   September 14, 2015, 04:24
Default
  #9
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
Hello,

I am working on similar project but for 3D wings.
The results of the ODE system are stored in y (in TestODE) which basically the initial conditions yStart.

Best wishes,
Hassan
meth likes this.
__________________
@HIKassem | HassanKassem.me
hk318i is offline   Reply With Quote

Old   October 7, 2015, 02:17
Default
  #10
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hi Hassan,

I went through the ODETest files and figured out how it works. I have implemented the ODE solver for the solid motion inside the icoFoam solver. Now, I need to take the force data to solve the solid motion solver in each time step. I went through the force.C and force.H files, but still have no idea how to use them to take the force data. I know the force.C file was written to write the force data, but I don't know how to edit those files to get the force data inside the icoFoam solver. I am not quite familiar with the library files and how OpenFoam uses the library files when we run a solver. So can you please help me on this?

Thanks.

Methma
meth is offline   Reply With Quote

Old   October 7, 2015, 07:17
Default
  #11
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
Quote:
Originally Posted by meth View Post
Hi Hassan,

I went through the ODETest files and figured out how it works. I have implemented the ODE solver for the solid motion inside the icoFoam solver. Now, I need to take the force data to solve the solid motion solver in each time step. I went through the force.C and force.H files, but still have no idea how to use them to take the force data. I know the force.C file was written to write the force data, but I don't know how to edit those files to get the force data inside the icoFoam solver. I am not quite familiar with the library files and how OpenFoam uses the library files when we run a solver. So can you please help me on this?

Thanks.

Methma
Hello Methma,

That is a great progress, well done.
Regarding the forces class, it is not only for writing. You can use it to get the forces directly [LINK].
There is a very good example for that in OpenFOAM source in sixDoFRigidBody [LINK] library. Here is a snippet of this library;

Code:
// define a dictionary first
    dictionary forcesDict;

// add the required parameters
    forcesDict.add("type", forces::typeName);
    forcesDict.add("patches", wordList(1, ptPatch.name()));
    forcesDict.add("rhoInf", rhoInf_);
    forcesDict.add("rhoName", rhoName_);
    forcesDict.add("CofR", motion_.centreOfRotation());

// create force object
    forces f("forces", db(), forcesDict);
// calculate the forces
    f.calcForcesMoment();
// finally you can get the forces and moment
    vector F = f.forceEff();
    vector M = f.momentEff();
Best wishes,
Hassan
meth likes this.
hk318i is offline   Reply With Quote

Old   October 7, 2015, 19:38
Default
  #12
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Thanks a lot Hassan, I will go through the links you mentioned. I hope this is very helpful.

Best,

Methma
meth is offline   Reply With Quote

Old   October 29, 2015, 01:22
Default
  #13
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hi Hassan,

I need your help again. I still could not do it. I have not quite understood how to deal with objectRegestry.

I know I have to give the corresponding objectRegestry as an argument to the force class when we construct a force object f. I have a little knowledge on objectRegistry class. I need to compute the force on the sphere boundary. So how should I initialise the objectRegestry and use it? Can you please give me some ideas on this.

I greatly appreciate your help.

Thank you.

Methma.
meth is offline   Reply With Quote

Old   October 29, 2015, 05:01
Default
  #14
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
Quote:
Originally Posted by meth View Post
Hi Hassan,

I need your help again. I still could not do it. I have not quite understood how to deal with objectRegestry.

I know I have to give the corresponding objectRegestry as an argument to the force class when we construct a force object f. I have a little knowledge on objectRegistry class. I need to compute the force on the sphere boundary. So how should I initialise the objectRegestry and use it? Can you please give me some ideas on this.

I greatly appreciate your help.

Thank you.

Methma.
You don't have to initialise an objectRegestry, you all need to pass the current objectRegestry to construct the forces object. The solution is quite easy, you can pass a reference to time or mesh [LINK]. Actually both are objectRegestry objects. Also, if you have access to any field like velocity, you can get reference to objectRegestry using db() function.

Bw,
Hassan
meth likes this.
hk318i is offline   Reply With Quote

Old   October 30, 2015, 03:37
Default
  #15
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
hi Hassan,

When I simple db() as an argument to the force object in the compilation of the solver it pop up the following error

vivicoFoam.C: In function ‘int main(int, char**)’:
vivicoFoam.C:137:27: error: ‘db’ was not declared in this scope
forces f("forces", db(), forcesDict);

It seems like it does not automatically take the db() function of the objectRegistry object. So, how can I fix that?

Thanks,

Methma
meth is offline   Reply With Quote

Old   October 30, 2015, 03:51
Default
  #16
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
If you have mesh or time reference pass it, if not try the velocity Field U.dB()
meth likes this.
__________________
@HIKassem | HassanKassem.me
hk318i is offline   Reply With Quote

Old   October 31, 2015, 23:47
Default
  #17
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Quote:
Originally Posted by hk318i View Post
If you have mesh or time reference pass it, if not try the velocity Field U.dB()

Thank you very much Hassan , It works, I think now things get clear and clear for me.

Best,

Methma
meth is offline   Reply With Quote

Old   November 1, 2015, 19:28
Default
  #18
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
Hassan,

Thank you very much, with your help I did it. Now I am able to calculate the forces on the solid. So, I can move into the next step of the solver.

Thanks again

Best,

Methma
hk318i likes this.
meth is offline   Reply With Quote

Old   November 2, 2015, 19:25
Default
  #19
Member
 
methma Rajamuni
Join Date: Jul 2015
Location: Victoria, Australia
Posts: 40
Rep Power: 10
meth is on a distinguished road
I want to impose a dirichlet boundary condition in Openfoam on the inlet patch with a time varying value calculated in each time step. For a example U = ( 1 y 0), where y is a value that I am calculating in each time step inside the solver. I know I have to use fixedValue type boundary condition which is uniform. Is there any inbuilt boundary condition which I can directly use or, do I need to create my own one? If I have to create a new boundary condition can you please give me some ideas and references.

Thank you

Methma
meth is offline   Reply With Quote

Old   November 5, 2015, 06:29
Default
  #20
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: Germany
Posts: 242
Rep Power: 17
hk318i is on a distinguished road
Hi Methma,

There are few options you have to tackle this case. First, you can create a new BC which calculates the velocity. For that you can check the time varying BCs. Your second option is to use codeStream to implement such BC. Your third option is to calculate this BC directly with the solver and updated each time step. This option discussed many times in this forum and you will find many useful tips on how to do it.
I will not go with the last option unless it is the only way because the BC will be solver dependent.

Bw,
Hassan
meth likes this.
hk318i is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


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


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