[ff3d-users] Variational form of potential flow with free surface on mesh

[Thomas Ward]

First of all, thanks to the ff3d team for making ff3d available, I found
the approach very natural and the documentation excellent.


I am using ff3d to solve a potential flow problem for forced periodic
oscillation of fluid in a tank with gravity free surface waves.
The solution using "strong" form appears to work, see below.


vertex a = (0,0,0);
vertex b = (10,10,10);
vector n = (25,50,10);

mesh M = structured(n,a,b);

double f=1.4;
function alpha=f*f/9.81;
solve(phi) in M
{
pde(phi)
-div(grad(phi))=0;
dnu(phi)=f on M xmin; //forced velocity = - omega cos(wt)
dnu(phi)=-f on M xmax;
dnu(phi)=0 on M ymin;// normal velocity zero on y walls
dnu(phi)=0 on M ymax;
dnu(phi)=0 on M zmin; //normal velocity zero on bottom
alpha*phi +dnu(phi) = 0 on M zmax; //free surface condition
};
save(medit,"tank1-phi",M);
save(medit,"tank-phi",phi,M);



I then tried to solve the same problem using the "weak" form

solve(phi) in M
{
test(w)
int (dx(w)*dx(phi)+dy(w)*dy(phi)+dz(w)*dz(phi))
// int(grad(phi)*grad(w))
+int[M xmin](w*dx(phi) +w*f)
+int[M xmax](w*dx(phi)-w*f)
+int[M ymin](w*dy(phi))
+int[M ymax](w*dy(phi))
+int[M zmin](w*dz(phi))
+int[M zmax](w*dz(phi)+w*alpha*phi)
=0;
};


but I get the following error

=========================
== Text mode execution ==
=========================
FreeFEM3D computing thread: started
Parsing data

VariationalFormulaExpression.cpp:73:
unexpected operator type

UNEXPECTED ERROR: this should not occure, please report it

BUG REPORT: Please send bug reports to:
ff3d-***@nongnu.org or ***@ann.jussieu.fr
or better, use the Bug Tracking System:
http://savannah.nongnu.org/bugs/?group=ff3d




any ideas?


PS

One small documentation suggestion; the method of defining Robin type
boundary conditions be added to the documentation or README file.


Tom

[Stephane Del Pino]

Dear Thomas,
thanks for your interest in ff3d.

In fact, you did not compute the right variational formulation associated to
your PDE. You have to be careful in the integration by part

Have a look to
http://en.wikipedia.org/wiki/Finite_element_method#Variational_formulation
and to
http://en.wikipedia.org/wiki/Green%27s_identities
for some details.

In your case, you should have written:

solve(phi) in M
{
test(w)
int(grad(phi)*grad(w))
+int[M zmax](w*alpha*phi)
=int[M xmin](w*f)
-int[M xmax](w*f)
=0;
};

Do not hesitate to ask me if you get troubles with that.

Post by Thomas Ward
PS
One small documentation suggestion; the method of defining Robin type
boundary conditions be added to the documentation or README file.

The documentation is far from being completed. I not this on the TODO list.

Best regards,
Stéphane.

[Thomas Ward]

Thanks Stephane, that fixed it.

I thought my formulation was equivalent to the one you wrote,

As I understand it, my formulation includes the natural boundary
conditions which is unnecessary but I thought that there was no harm
(other than inefficiency) to include them, anyway thanks for the
pointers to the wikipedia stuff, I will look through them and try to see
why including the natural boundary codtions in the surface integral
won't work.


I'm about to try using periodic BC, I note from the mailing list that
you were playing around with this a couple of years ago, do you have any
sample files or other documentation or pointers to the source code?

Do you have any infinite elements? I'm thinking of a radiation condition
decaying to zero at infinity on an open boundary.

regards,

Tom

Post by Stephane Del Pino
Dear Thomas,
thanks for your interest in ff3d.
In fact, you did not compute the right variational formulation associated to
your PDE. You have to be careful in the integration by part
Have a look to
http://en.wikipedia.org/wiki/Finite_element_method#Variational_formulation
and to
http://en.wikipedia.org/wiki/Green%27s_identities
for some details.
solve(phi) in M
{
test(w)
int(grad(phi)*grad(w))
+int[M zmax](w*alpha*phi)
=int[M xmin](w*f)
-int[M xmax](w*f)
=0;
};
Do not hesitate to ask me if you get troubles with that.

Post by Thomas Ward
PS
One small documentation suggestion; the method of defining Robin type
boundary conditions be added to the documentation or README file.

The documentation is far from being completed. I not this on the TODO list.
Best regards,
Stéphane.
_______________________________________________
ff3d-users mailing list
http://lists.nongnu.org/mailman/listinfo/ff3d-users

[Stephane Del Pino]

Hello Thomas.

Post by Thomas Ward
Thanks Stephane, that fixed it.
I thought my formulation was equivalent to the one you wrote,
As I understand it, my formulation includes the natural boundary
conditions which is unnecessary but I thought that there was no harm
(other than inefficiency) to include them, anyway thanks for the
pointers to the wikipedia stuff, I will look through them and try to see
why including the natural boundary codtions in the surface integral
won't work.

It is not equivalent. I join you a formal way of getting this variationnal
formula. Check finite element books to get details and a clean establishment
of the formula.

Post by Thomas Ward
I'm about to try using periodic BC, I note from the mailing list that
you were playing around with this a couple of years ago, do you have any
sample files or other documentation or pointers to the source code?

Yes this is not a big deal in ff3d. Assume that M is your structured mesh,
then
M = periodic(M,0:1,2:3,4:5);
will create a triperiodic mesh where
0 represents xmin,
1 represents xmax,
2 represents ymin,
3 represents ymax,
4 represents zmin,
5 represents zmax,

M = periodic(M,2:3);
would create a periodic mesh in the direction y. Note that up to now M must be
a structured mesh.

Post by Thomas Ward
Do you have any infinite elements? I'm thinking of a radiation condition
decaying to zero at infinity on an open boundary.

No. No such element is implemented ...

Best regards,
Stéphane.

[Thomas Ward]

Many thanks Stephane, I think I have got it now.

we choose test_fn_boundary = -test_fn_volume so that after integrating
the volume integral of test_fn_volume x div(grad(phi) by parts the
surface integral terms cancel out with Neumann terms of the
test_fn_boundary x boundary conditions integral.


I hope to try the periodic conditions next week, its not totally clear
to me how the periodic mesh becomes a periodic boundary condition but I
will have a try first and then ask if I can't work it out.


Tom

Post by Stephane Del Pino
Hello Thomas.

Post by Thomas Ward
Thanks Stephane, that fixed it.
I thought my formulation was equivalent to the one you wrote,
As I understand it, my formulation includes the natural boundary
conditions which is unnecessary but I thought that there was no harm
(other than inefficiency) to include them, anyway thanks for the
pointers to the wikipedia stuff, I will look through them and try to see
why including the natural boundary codtions in the surface integral
won't work.

It is not equivalent. I join you a formal way of getting this variationnal
formula. Check finite element books to get details and a clean establishment
of the formula.

Post by Thomas Ward
I'm about to try using periodic BC, I note from the mailing list that
you were playing around with this a couple of years ago, do you have any
sample files or other documentation or pointers to the source code?

Yes this is not a big deal in ff3d. Assume that M is your structured mesh,
then
M = periodic(M,0:1,2:3,4:5);
will create a triperiodic mesh where
0 represents xmin,
1 represents xmax,
2 represents ymin,
3 represents ymax,
4 represents zmin,
5 represents zmax,
M = periodic(M,2:3);
would create a periodic mesh in the direction y. Note that up to now M must be
a structured mesh.

Post by Thomas Ward
Do you have any infinite elements? I'm thinking of a radiation condition
decaying to zero at infinity on an open boundary.

No. No such element is implemented ...
Best regards,
Stéphane.
------------------------------------------------------------------------
_______________________________________________
ff3d-users mailing list
http://lists.nongnu.org/mailman/listinfo/ff3d-users

[Thomas Ward]

Stephane

last May you gave me a pointer to how to work with periodic boundary
conditions in ff3d. 8 months later I have found the time to try it.....

1.
I'm not able to get it to work, my file is attached at the end of this
message can you see what I am doing wrong?

2.
I would like to understand how it is implemented in the code have been
trying to understand code but can't manage to follow it. Can you point
me in the right directions, so far I have been looking through
MeshPeriodizer.cpp and VerticesCorrespondance.cpp , Mesh.cpp and
SurfaceMeshOfQuadrangles.cpp

I can see that the nodes of the end boundary of the structured mesh are
mapped to the nodes of the start boundary but I am struggling to see how
this is carried forwards to the solver.

Are there any papers which describe the method you have followed?



thanks for your help

Tom Ward.



vertex a = (0,0,0);
vertex b = (10,10,10);
vector n = (25,50,10);

mesh M = structured(n,a,b);
//this next line is the only change I made to the static solution
M = periodic(M,0:1);

double f=1.4;
function alpha=f*f/9.81;
solve(phi) in M
{
test(w)
int(grad(phi)*grad(w))
+int[M zmax](w*alpha*phi)
+int[M xmin](w*f)
-int[M xmax](w*f)
=0;
};

save(vtk,"tank-varA-periodic",phi,M);

Post by Thomas Ward
Many thanks Stephane,
I hope to try the periodic conditions next week, its not totally clear
to me how the periodic mesh becomes a periodic boundary condition but I
will have a try first and then ask if I can't work it out.
Tom

......

Post by Thomas Ward

Post by Stephane Del Pino

Post by Thomas Ward
I'm about to try using periodic BC, I note from the mailing list that
you were playing around with this a couple of years ago, do you have any
sample files or other documentation or pointers to the source code?

Yes this is not a big deal in ff3d. Assume that M is your structured
mesh, then
M = periodic(M,0:1,2:3,4:5);
will create a triperiodic mesh where
0 represents xmin,
1 represents xmax,
2 represents ymin,
3 represents ymax,
4 represents zmin,
5 represents zmax,
M = periodic(M,2:3);
would create a periodic mesh in the direction y. Note that up to now M
must be a structured mesh.

[Stephane Del Pino]

Dear Thomas.

I think that the problem comes from the fact that you impose boundary
conditions on xmin and xmax and also mark these boundary as being periodic.
I just checked and periodic boundary conditions still work.

Best regards,
Stéphane.