[Gfs-devel] numerical integration

[Daniel Fuster]

Hi Stephane

As you know, I would like to extend the capabilities of
http://gfs.sourceforge.net/wiki/index.php/GfsOutputErrorNorm
in order to allow for more integration points when we integrate some
analytical function on the cell

You told me once that you had coded already some functions that may be
useful for this end. I would be very grateful if you could share them.
The other possibility is to implement Gaussian quadrature and to
specify the coefficients for low order quadratures inside the code. It
should not be too hard to do this

PS: I tried to look if there were something in GSL but I only found
functions for 1D integration...

Thanks
Daniel

------------------------------------------------------------------------------
Everyone hates slow websites. So do we.
Make your web apps faster with AppDynamics
Download AppDynamics Lite for free today:
http://p.sf.net/sfu/appdyn_d2d_nov
_______________________________________________
Gfs-devel mailing list
Gfs-...@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/gfs-devel

[Geordie McBain]

2012/11/10 Daniel Fuster <dfu...@gmail.com>:

> Hi Stephane
>
> As you know, I would like to extend the capabilities of
> http://gfs.sourceforge.net/wiki/index.php/GfsOutputErrorNorm
> in order to allow for more integration points when we integrate some
> analytical function on the cell
>
> You told me once that you had coded already some functions that may be
> useful for this end. I would be very grateful if you could share them.
> The other possibility is to implement Gaussian quadrature and to
> specify the coefficients for low order quadratures inside the code. It
> should not be too hard to do this
>
> PS: I tried to look if there were something in GSL but I only found
> functions for 1D integration...

Hi Daniel. I'm guessing that the omission from GSL might be because
the most common thing to do on rectangles is use what's called a
"tensor-product" of the one-dimensional rule. I think this how most
finite element codes evaluate inner products on their "quads" and
"hexes". Basically if int_0^1 f (x) dx is approximated by sum_i
f(x_i) w_i then int_0^1 int_0^1 f (x, y) dx dy is approximated by
sum_{i,j} w_i f (x_i, x_j) w_j. There are specialized
multidimensional rules, but they're somewhat arcane; see Table 6.7.2
of Engels, H. (1980, Numerical Quadrature and Cubature, Academic) for
a third-order three-point formula for the square.

But perhaps another approach making more use of existing Gerris
infrastructure would be an "adaptive quadrature" (see
http://en.wikipedia.org/wiki/Adaptive_quadrature) based on octrees.
Although not as theoretically optimal as Gaussian quadrature in
one-dimension (and tensor product Gaussian quadrature isn't optimal in
two or more dimensions anyway), adaptive quadrature tends to be more
robust in cases of sharp localized variation. A test which should
highlight that would be to compute the volume of a circle or sphere by
integrating its characteristic function (f(x) = inside ? 1 : 0),
though perhaps that's not a fair test as you did specify that the
integrand be analytic.
I guess the code implementing that would be largely a combination of
GfsOutputScalarSum (or GfsOutputErrorNorm) and GfsAdapt.

[Daniel Fuster]

Thanks a lot for your comments Geordie

I think we will give a try to "adaptive quadrature". As you say, the
only hard part is that it will be required to create new virtual
"refined" cells to obtain accurate values of the integral. The
difference with adapt is that it will be necessary to create cells
that may be several levels of refine "deeper" that the current leaf
cell (and also destroy them once the integral has been obtained).

The idea is to include new optional arguments in OutputErrorNorm, namely
maxlevel : which can be larger than the maxlevel in domain
minlevel: the minimum level in which the integration will be done
tolerance: Stop criterion to adapt the mesh if the value of the
integral has converged from one level to next level

Thanks
Daniel

2012/11/12 Geordie McBain <gdmc...@freeshell.org>:

--
www.danielfuster.com

------------------------------------------------------------------------------
Monitor your physical, virtual and cloud infrastructure from a single
web console. Get in-depth insight into apps, servers, databases, vmware,
SAP, cloud infrastructure, etc. Download 30-day Free Trial.
Pricing starts from $795 for 25 servers or applications!
http://p.sf.net/sfu/zoho_dev2dev_nov

[Daniel Fuster]

Hi

I think I have a preliminary patch that do something useful. Before
going into adaptive integration, I have done a patch that allows to
integrate the error at an arbitrary level of refinement (for now at
constant level). This will allow for accurate error norm calculations
if desired.

In the attached patch, I create new childs when required in order to
do the integration at the required level. Then the norm is computed
and finally the new cells are removed. I had to add some general
functions in the code that I hope are self-explained.

The patch should give exactly the same answer than the previous
OutputErrorNorm in cases where maxlevel is not specified (which I
think is always the case). When the user specifies its value, the
integration is done at the user-defined level. I attach one example
where the error with respect to the reference solution is the area of
a circle. In this case the exact solution is pi/100. One gets the
right answer with a level=1 simulation doing:

gerris2D -m -DLEVEL=1 -DMAXLEVEL=10 integral.gfs

Eventually one could do also adaptive integration with some
modifications. For simplicity I prefer to send this patch and wait for
comments.

PS: The error calculation applying AMR will be different with respect
to the stable release if one specifies the finest level in
outputerrornorm because all the solution will be projected in a fine
grid to obtain the error

Best
Daniel

1 0 GfsSimulation GfsBox GfsGEdge {} {
Global {
}

Time { iend = 0 dtmax = 1.e-3 }

GfsVariableTracerVOF T

GfsInit {} { U = 0. }

Refine LEVEL

OutputErrorNorm { start = end } stdout {
v = U
maxlevel = MAXLEVEL
} {
# s = exp(-(x*x+y*y)/0.005)
s = (x*x + y*y - 0.1*0.1 < 0. )
v = Error
}

}
GfsBox {}


2012/11/12 Daniel Fuster <dfu...@gmail.com>:

www.danielfuster.com

[Stephane Popinet]

Hi Daniel,

Could you please explain a bit more the background rationale for
higher-order integration in OutputErrorNorm? I am not sure I
understand why this is required.

The routines I was mentioning use "sparse grids" and are in
metric.c:209 together with links to reference papers.

cheers

Stephane

------------------------------------------------------------------------------
Monitor your physical, virtual and cloud infrastructure from a single
web console. Get in-depth insight into apps, servers, databases, vmware,
SAP, cloud infrastructure, etc. Download 30-day Free Trial.
Pricing starts from $795 for 25 servers or applications!
http://p.sf.net/sfu/zoho_dev2dev_nov

[Daniel Fuster]

Hi Stephane

>
> Could you please explain a bit more the background rationale for
> higher-order integration in OutputErrorNorm?

Some error estimation techniques are based on the exact integration of
the error function, defined as the difference between the exact and
numerical solution. The error function can be a high order polynomial
for the general case that requires high order integration techniques.

Generally, high order integration can provides (sometimes) more
accurate and also consistent measures of the error. This is especially
true when looking at local errors. The reason is that if one applies
the mid point rule at the cell center, one can see small errors by
chance, that misbehave when the mesh is refined

For instance, in the simple example I included, if I assume that the
solution provided by the numerical method with small resolution is 0
(say at level 1), and we want to compute the error with respect the
real solution (assumed to be 1 inside the circle and 0 outside), one
gets a very small error (even 0!) if the integration is done at level
1.

Just to make things clear, this is only an alternative way of
computing the error (not universal) that it is sometimes useful to
consider in order to investigate the quality of the solution. I think
that allowing the user to choose which kind of integration wants to
perform is useful.

Please, tells me if that makes sense to you
Daniel

>I am not sure I
> understand why this is required.
>
> The routines I was mentioning use "sparse grids" and are in
> metric.c:209 together with links to reference papers.
>
> cheers
>
> Stephane
>
> ------------------------------------------------------------------------------
> Monitor your physical, virtual and cloud infrastructure from a single
> web console. Get in-depth insight into apps, servers, databases, vmware,
> SAP, cloud infrastructure, etc. Download 30-day Free Trial.
> Pricing starts from $795 for 25 servers or applications!
> http://p.sf.net/sfu/zoho_dev2dev_nov
> _______________________________________________
> Gfs-devel mailing list
> Gfs-...@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/gfs-devel

--
www.danielfuster.com

[Stephane Popinet]

Hi Daniel,

> Please, tells me if that makes sense to you

I think I see where you are coming from, however I am not sure about
the details. For example, evaluating the integral over a cell to
higher-order involves the knowledge of the numerical solution at
arbitrary points in the cells (or equivalently the coefficients of
basis functions defining a continuous solution over the entire cell).
This interpolation should be consistent with the underlying scheme.
There is currently no obvious way to do this consistently in the code.

[Vladimir Kolobov]

Daniel,

are you talking about calculating integrals over arbitrary (predefined)
surfaces for a function defined at cell centers ?

Thanks,

Vladimir

[Daniel Fuster]

> This interpolation should be consistent with the underlying scheme.

You are enterely right. For now, I consider the numerical solution to
be constant within the cell, which is a general assumption for error
estimation techniques applied to finite volumes. As you say, this is
not rigorously true for gerris because it uses "special" interpolation
proceduce to obtain, for example, the velocity at the faces. Indeed,
in gerris, one could say that the representation of the numerical
solution depends on the variable considered (the numerical scheme is
not the same for a tracer, than for vof, than for momentum...).

Despite of these uncertitudes, I still think this approach can give in
some cases a better approximation of the error than the mid point
rule applied to the leaf cells (better approximation does not mean
exact calculation). In a general case, it will be useful sometimes and
some other times it will make no sense to use it. To allow the user to
have the possibility to use it I think it is possitive thing. Taking
my example only as a particular case, the method I propose gives a
better measure of the "quality" of the representation of the numerical
solution with respect to the real solution at low resolution levels.
This happens for this particular case (and also with many others) but
I am sure that you can find plenty of counter-examples. Again, giving
the user the choice of getting two different "approximate" measures of
the error I think it is a good thing

Vladimir, in principle I just want to compute the error of a given
function with respect to a cell centered solution over the whole
domain. For now I do not consider integration over pre-defined
surfaces (but maybe the discussion applies also)

best
Daniel

> There is currently no obvious way to do this consistently in the code.
>
> cheers
>
> Stephane
>
> ------------------------------------------------------------------------------
> Monitor your physical, virtual and cloud infrastructure from a single
> web console. Get in-depth insight into apps, servers, databases, vmware,
> SAP, cloud infrastructure, etc. Download 30-day Free Trial.
> Pricing starts from $795 for 25 servers or applications!
> http://p.sf.net/sfu/zoho_dev2dev_nov
> _______________________________________________
> Gfs-devel mailing list
> Gfs-...@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/gfs-devel

--
www.danielfuster.com

[Geordie McBain]

2012/11/13 Vladimir Kolobov <v...@cfdrc.com>:

> Daniel,
>
> are you talking about calculating integrals over arbitrary (predefined)
> surfaces for a function defined at cell centers ?

No. That sounds interesting, but current interest is restricted to the
domain. The context is a posteriori error estimation, and in
particular assessing schemes for that using test cases with known
analytic solutions to see how efficient the estimator is in each cell.

[Daniel Fuster]

Hi Stephane

I just send another patch that does not require any patch related to
check_neighbors flags.

Also an example where I try to get a measure of the accuracy of the
representation of a VOF field for different levels of refinement. For
this example, the high order integration gives a more consistent
representation of the "accuracy" of the solution at low levels of
refinement. To run the example just type

sh errorcircle.sh

To have various and good error measurements (also locally) is very
useful when one wants to optimize AMR

PS: Thanks for your effort of replying all the pending emails

Daniel

2012/11/13 Daniel Fuster <dfu...@gmail.com>:

www.danielfuster.com

[Stephane Zaleski]

Hi Daniel

I have some trouble understanding what you do, and that reflects our discussion
with Geordie in my office and the things I have written on the whiteboard.

First of all I think that you should define the error function as the
"residual" and not the "difference between the exact and
numerical solution".

Second are you dealing with these cases where the exact solution is
somehow known, for example analytically ?

Third, in the VOF case the solution is not differentiable, actually not even
continuous, which might affect things a great deal ?

Finally, exact initialization of VOF functions is something Ruben Scardovelli
is currently working on so perhaps it would be interesting to email him ?
he will be in Paris in January.

best

Stéphane

PS I wrote to Geordie to ask him formally what he was working on until his resignation
takes effect and he told me he was working with you on these error estimation things.
I hope you find the help useful.

For our TOTAL project we desperately need some improvement in the speed of Gerris
and probably a better way to represent solids too. I hope that bertrand is interacting with you.
I would gladly help but I have my hands full with ParisSimulator at the moment. When that
has made progress I will be able to help with Gerris again.

[Stephane Zaleski]

I am sorry this message was not meant for the mailing list but for Daniel personnaly
please disregard it.