Integration(math) issue

sabago · November 9, 2015, 15:50

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

alexeym · November 10, 2015, 03:20

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.

sabago · November 10, 2015, 03:51

Quote:

Originally Posted by alexeym

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?

alexeym · November 10, 2015, 04:10

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

sabago · November 10, 2015, 04:27

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"

alexeym · November 10, 2015, 04:37

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".

sabago · November 10, 2015, 04:41

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

alexeym · November 10, 2015, 04:58

OK, your case is 2D, so is jc a surface? So you have

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?

sabago · November 10, 2015, 05:03

Quote:

Originally Posted by alexeym

OK, your case is 2D, so is jc a surface? So you have

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

alexeym · November 10, 2015, 05:30

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.

November 10, 2015, 03:20		#2
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
[swak4Foam] swak4Foam 0.3.1 OF2.2.2 Installation Issue	scnsc	OpenFOAM Community Contributions	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 02: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

November 10, 2015, 04:27		#5
sabago Member Sandra Join Date: Oct 2014 Posts: 58 Rep Power: 11	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"

November 10, 2015, 04:37		#6
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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".

November 10, 2015, 04:41		#7
sabago Member Sandra Join Date: Oct 2014 Posts: 58 Rep Power: 11	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

November 10, 2015, 04:58		#8
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	OK, your case is 2D, so is jc a surface? So you have 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?

November 10, 2015, 05:30		#10
alexeym Senior Member Alexey Matveichev Join Date: Aug 2011 Location: Nancy, France Posts: 1,930 Rep Power: 38	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 ax + b. You can use polyLineSet set class to construct cell set. 3. You mesh moves, is irregular, your line is defined as ax + b. You reconstruct cell set using polyLineSet after every mesh topological change.