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

Integration(math) issue

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

Reply
 
LinkBack Thread Tools Display Modes
Old   November 9, 2015, 16:50
Default Integration(math) issue
  #1
Member
 
Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
sabago is on a distinguished road
Hello Foamers!
This code might look familiar but I have added a few lines that require an integration that isn't going so well. Anyhoo, here's the code

Code:
forAll(x3, celli)
{
for (count=0; count<iter; count++) 
{
        x3 = (x1 + x2)/2;

       jc = i*(Foam::exp(8.1*x))-(Foam::exp(-2.7*x))

	f1 = (Foam::exp(8.1*x1))-(Foam::exp(-2.7*x1))-(jc/i);

	f3 = (Foam::exp(8.1*x3))-(Foam::exp(-2.7*x3))-(jc/i);
	
	 if(neg(f1[celli] *f3[celli]))
        {
            x2 = x3;
        }
        else if((mag(x1[celli] - x2[celli]) < 1e-7)||(f3[celli]==0)||(jc[celli].integrate(0,3.6e9)==2e4))
        {
            root =x3;
	    break;
        }
	else
	{
	x1= x3;
	}
}
     
Info<<"Root is"<<endl;
Info<<root<<endl;
}
jc is declared as a volScalarField in my createFields and here's the error that I get (I understand that I do not have a double ie I understand what the error means but I do not know how to proceed.

Quote:
error: request for member ‘integrate’ in ‘jc.Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::<anonymous>.Foam:imensionedField <double, Foam::volMesh>::<anonymous>.Foam::Field<double>::< anonymous>.Foam::List<double>::<anonymous>.Foam::U List<T>:perator[]<double>(celli)’, which is of non-class type ‘double’
Many thanks in advance,
Sandra
sabago is offline   Reply With Quote

Old   November 10, 2015, 04:20
Default
  #2
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,438
Rep Power: 25
alexeym will become famous soon enoughalexeym will become famous soon enough
Hi,

1.
Code:
f3[celli]==0
in general it is not a good idea to compare doubles with ==. There are two functions:

Code:
inline bool equal(const Scalar& s1, const Scalar& s2)
{
    return mag(s1 - s2) <= ScalarVSMALL;
}

    
inline bool notEqual(const Scalar s1, const Scalar s2)
{
    return mag(s1 - s2) > ScalarVSMALL;
}
2. if jc is volScalarField then jc[celli] is just scalar (i.e. simple number), so

\int_a^b Cdx = Cx|_a^b = C(b -a)



And again this

Code:
jc[celli].integrate(0,3.6e9)==2e4
will be FALSE 99,9% of the time.
alexeym is offline   Reply With Quote

Old   November 10, 2015, 04:51
Default
  #3
Member
 
Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
sabago is on a distinguished road
Quote:
Originally Posted by alexeym View Post
2. if jc is volScalarField then jc[celli] is just scalar (i.e. simple number), so

\int_a^b Cdx = Cx|_a^b = C(b -a)



And again this

Code:
jc[celli].integrate(0,3.6e9)==2e4
will be FALSE 99,9% of the time.
Would that formula work even if C is not a linear function per say (and I don't know its form either)? There's no method in OpenFOAM to get the area under a curve during simulation for volScalarFields like there's for doubles ie the integrate that I tried to use?

Lastly, out of curiosity, why would
Code:
jc[celli].integrate(0,3.6e9)==2e4
be false 99.9% of the time?
sabago is offline   Reply With Quote

Old   November 10, 2015, 05:10
Default
  #4
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,438
Rep Power: 25
alexeym will become famous soon enoughalexeym will become famous soon enough
If, as you said, jc is volScalarField, then jc[celli] can not be "...not a linear function per say...". It is a number and on a given time step it is in general even constant. I guess, you need to explain what you are trying to achieve with jc[celli].integrate(...). Right now what you are trying to do is like this:

Code:
(3.1415).integrate(0, 3e9)
Why do you need to integrate constant? What is integration variable?

Quote:
Lastly, out of curiosity, why would
...
be false 99.9% of the time?
https://isocpp.org/wiki/faq/newbie#floating-point-arith
https://en.wikipedia.org/wiki/Floati...uracy_problems
https://duckduckgo.com/?q=floating+point+comparison
alexeym is offline   Reply With Quote

Old   November 10, 2015, 05:27
Default
  #5
Member
 
Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
sabago is on a distinguished road
sorry, my bad....

jc in paraview (at each time step) is a curve like a straight up witch hat and I want to iterate until the area under that curve equals a constant that I give it (in this case, 20,000).
So, I agree, jc[celli] wouldn't work here but I wasn't thinking of that, unfortunately, when I tried it after "jc.integrate" gave me an error along the lines of "member doesn't have function integrate"
sabago is offline   Reply With Quote

Old   November 10, 2015, 05:37
Default
  #6
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,438
Rep Power: 25
alexeym will become famous soon enoughalexeym will become famous soon enough
Well, let us start from the beginning...

What is jc? In the first post you wrote, that it is volScalarField, i.e. something that varies in VOLUME, then you have written, that it is "...a curve like a straight up witch hat...", i.e. someting that varies just in ONE direction (surely there are 3D curves, yet they usually can be reduced to function of one variable).

After you define jc, I think, it is easier to invent a method for calculation your "integral".
alexeym is offline   Reply With Quote

Old   November 10, 2015, 05:41
Default
  #7
Member
 
Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
sabago is on a distinguished road
it is defined as a volScalarField and calculated from the equation in the initial post BUT my case is 2D ie I have front and back patches empty so technically it's an areaScalarField if that were a thing in OpenFOAM and that's why I talked about a 2d curve.
Does that help?

Sandra
sabago is offline   Reply With Quote

Old   November 10, 2015, 05:58
Default
  #8
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,438
Rep Power: 25
alexeym will become famous soon enoughalexeym will become famous soon enough
OK, your case is 2D, so is jc a surface? So you have

jc = jc(x, y)

and would like to calculate

\int_{x_a}^{x_b} \int_{y_a}^{y_b} jc\left(x, y\right)dxdy

where xa, xb, ya, yb are the corners of your simulation domain (so it would be volume under the curve and in simplest form it is a sum of jc[celli]*mesh.V()[celli]/<distance in empty direction>)?

Or you would like to take values of jc along the line L (technically reducing surface to curve) and them integrate 1D function along L?
alexeym is offline   Reply With Quote

Old   November 10, 2015, 06:03
Default
  #9
Member
 
Sandra
Join Date: Oct 2014
Posts: 58
Rep Power: 4
sabago is on a distinguished road
Quote:
Originally Posted by alexeym View Post
OK, your case is 2D, so is jc a surface? So you have

jc = jc(x, y)

and would like to calculate

\int_{x_a}^{x_b} \int_{y_a}^{y_b} jc\left(x, y\right)dxdy

where xa, xb, ya, yb are the corners of your simulation domain (so it would be volume under the curve and in simplest form it is a sum of jc[celli]*mesh.V()[celli]/<distance in empty direction>)?

Or you would like to take values of jc along the line L (technically reducing surface to curve) and them integrate 1D function along L?
The second option, please
\int_{x_a}^{x_b}jc\left(x \right)dx

Many thanks, BTW
sabago is offline   Reply With Quote

Old   November 10, 2015, 06:30
Default
  #10
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,438
Rep Power: 25
alexeym will become famous soon enoughalexeym will become famous soon enough
There are several possibilities:

1. Your mesh is static regular orthogonal, your line is parallel to axis. You construct cell set (cells, where you line passes) yourself in the beginning of the simulation, starting from xa (using, for example, findCell method of mesh) and then iterating through the neighbors and selecting cells in required direction (or if mesh is not large, you can use findCell for every point of your line, as you are doing it just once). After in simplest case you iterate through the labels, take value of jc in the cell, take width of the cell and multiply.

2. You mesh is static irregular, your line is defined as a*x + b. You can use polyLineSet set class to construct cell set.

3. You mesh moves, is irregular, your line is defined as a*x + b. You reconstruct cell set using polyLineSet after every mesh topological change.
alexeym 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
swak4Foam 0.3.1 OF2.2.2 Installation Issue scnsc OpenFOAM Installation 3 August 20, 2014 07:59
Strange issue while launching ANSYS workbench in CentOS 6.4 Philip_C ANSYS 11 August 29, 2013 06:44
CyclicAMI Issue In OpenFOAM 2.2.0 prasant OpenFOAM Running, Solving & CFD 17 March 16, 2013 03:00
Pressure boundary condition issue Vijay FLUENT 0 April 6, 2012 13:35
Meshing related issue in Flow EFD appu FloEFD, FloWorks & FloTHERM 1 May 22, 2011 08:27


All times are GMT -4. The time now is 00:42.