embed solid with metric + vof + solvers

[jose.lopez...@gmail.com]

Hi Stephane,

I have  advanced in my goal of making compatible metric + vof + embed + solvers.

 

For the moment:

1.- As you suggested fms now gathers either metrics and embed factors.

2.- VOF can interact with solids preserving the volume (no tricks would be necessary).

3.- Poisson solver now is compatible with metrics (for the moment just axi.h)

Now I want to adapt the viscous solver to deal either with embed and metric. 

 As you know, the version you use ("viscosity-embed.h") neglects some viscous terms that would appears if mu were not constant (as it would be in a two phase problem) . On the contrary "viscosity.h" keeps all terms.

I am not sure  what would be the best, to keep on viscosity-embed.h making the proper modifications for metric  despite some terms are lost, or center the efforts in "viscosity.h" adding the embed stuff. What do you think?

Also, could you explain me why you do:

p.embed_flux = u.x.boundary[embed] != antisymmetry ? embed_flux : NULL;

Cheers

Pepe

[Stephane Popinet]

Hi Pepe,

> 1.- As you suggested fms now gathers either metrics and embed factors.
> 2.- VOF can interact with solids preserving the volume (no tricks would
> be necessary).
> 3.- Poisson solver now is compatible with metrics (for the moment just
> axi.h)

That sounds excellent. Can you share the code and I will have a look?

I had started looking into this with the goal of making the
examples/tangaroa.c work with embedded boundaries. One issue was the
compatibility of the treatment of embed + VOF when combined with the
navier-stokes/conserving.h option, which is necessary for stability of
this air/water configuration. I had a first rough version which worked.
Do you think your version would be compatible with
navier-stokes/conserving.h?

> I am not sure  what would be the best, to keep on viscosity-embed.h
> making the proper modifications for metric  despite some terms are lost,
> or center the efforts in "viscosity.h" adding the embed stuff. What do
> you think?

I think starting from viscosity.h is probably the best option, but your
vision of this is probably "fresher" than mine!

> Also, could you explain me why you do:
>
> p.embed_flux = u.x.boundary[embed] != antisymmetry ? embed_flux : NULL;

Hmm, one important thing to do for viscosity*.h would be to write
documentation!

I am guessing that this has to do with whether default boundary
conditions (i.e. "symmetry") are set for u.x on the boundary. If this is
the case, then this corresponds to pure slip (i.e. Neumann) on the
embedded boundary and the "viscous" flux through the boundary is zero,
otherwise Dirichlet conditions are imposed and the flux must be computed.

cheers,

Stephane

[jose.lopez...@gmail.com]

Hi Stephane,

El miércoles, 24 de febrero de 2021 a las 9:12:37 UTC+1, Stephane Popinet escribió:

Hi Pepe,

That sounds excellent. Can you share the code and I will have a look?

 

Attached!

Do you think your version would be compatible with
navier-stokes/conserving.h?

 

I have not though in that issue since now I am focused in the viscous solver but I do not see any reason it should be incompatible.

> I am not sure  what would be the best, to keep on viscosity-embed.h
> making the proper modifications for metric  despite some terms are lost,
> or center the efforts in "viscosity.h" adding the embed stuff. What do
> you think?

I think starting from viscosity.h is probably the best option, but your
vision of this is probably "fresher" than mine!

 

Ok, annoted! I asked you because I wondered if you had some impelling reason (beyond a larger simplicity) to go back to a simpler viscous solver instead to the modern/complete version you have in viscosity.h

> Also, could you explain me why you do:
>
> p.embed_flux = u.x.boundary[embed] != antisymmetry ? embed_flux : NULL;

Hmm, one important thing to do for viscosity*.h would be to write
documentation!

Your are right I will document the files

Thanks, Stephane!

The masked  Pepe

[jmlopez...@gmail.com]

Hi  Stephane,

I forgot to tell you that  my patch produces a silly (and annoying) warning  that I have not been able to remove

Pepe

/src/embed.h:844:2: warning: control reaches end of non-void function [-Wreturn-type]

[Stephane Popinet]

Hi Pepe,

I had a quick look at the patch, it looks good and I will check it more
thoroughly later. A quick question though: it seems that many of the
changes are just cm -> cms. I can understand why this is so but did you
check which combination minimised the number of such changes? What I
mean is that if you swapped cm and cms (which you could then rename to
something more appropriate), would that lead to much less changes in the
code?

cheers,

Stephane

[Karl Cardin]

This in likely a naive question, but in the neumann-axi.c test case why does the right hand side change even though the exact solution in the the same as the Cartesian test? 

[Jose M. López-Herrera Sánchez]

in case someoneelse  want an answer...

---------- Forwarded message ---------
De: Jose M. López-Herrera Sánchez <jose.lopez...@gmail.com>
Date: mar, 2 mar 2021 a las 23:27
Subject: Re: embed solid with metric + vof + solvers
To: Karl Cardin <nidra...@gmail.com>

Hi Karl,

The laplace equation  is not the same in cylindrical coordinates and cartesian coordinates. In order to have the same solution necessarily the right hand side (b) should not be the same.

Cheers

Jose

--
You received this message because you are subscribed to a topic in the Google Groups "basilisk-fr" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/basilisk-fr/ul8K6HJhHzo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to basilisk-fr...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/basilisk-fr/1e202c71-9de1-4c90-b3c1-0c3eba5a6e2en%40googlegroups.com.

[Karl Cardin]

Hi Jose,

Thank you for your response. You cleared up my confusion.

And thanks for your work on metric + vof + solvers.

Best,

Karl

[Jose M. López-Herrera Sánchez]

Hi Stephane,

Finally, I cleaned the patch  and kept the changes to minimal (It took longer that I thought!). I have constructed the patch over the last release of basilisk,  I hope it serves to you.

Cheers

Pepe

PS. It is not very relevant but in my new workstation the test "hydrostatic2.c" fails writing in the stderr terminal.

El jue, 25 feb 2021 a las 15:57, Stephane Popinet (<pop...@basilisk.fr>) escribió:

> Shall I prepare my patch with (cm, fm) instead of (cms, fms)?

Yes, please.

Stephane

[xiaowei]

Dear Prof. Jose,

Thank you for the great contribution to the embed metric + vof.

After applying your patch, I did a few simple tests of a drop splashes on a wall with the wall defined by the embed function. It works fine with the following settings:

VOF + embed+ quadtree + MPI

VOF + embed + multigrid2D + MPI

VOF + embed + multigrid3D + MPI

VOF + embed + octree + MPI

VOF + embed + octree + (No) MPI

However, the VOF + embed + octree + MPI encountered divergence issues at the moment when the drop touched the embedded wall. The error report is as follows: 

WARNING: convergence for p not reached after 100 iterations

  res: 1.50561 sum: 5598.91 nrelax: 100

injection: /home/xiaozhu/basilisk/src/timestep.h:10: timestep: Assertion `((double *) ((((Tree *)grid)->L[point.level]->m[point.i+0][point.j+0]) + sizeof(Cell)))[(fm.x.i)]' failed.

So I’m wondering are there any compatible issues here, and is it possible to fix the issue? 

Best,

Xiaowei

[H LIU]

Hi Xiaowei,

Did your tests include "tension.h"? 

I test the patch too, but failed, and my test : embed+vof+surface tension.

Thanks.

Hu LIU

[ed mumei]

Dear José,

Thank you very much for this patch.

I tried to use it with the capillary wave test. I added the following lines to capwave.c:

(...)

#include "embed.h"

(...)

u.n[embed] = dirichlet(0.);
u.t[embed] = dirichlet(0.);
(...)

  vertex scalar phi[];
  foreach_vertex() {
    phi[] = sq(x-1.) + sq(y-0.5) - sq(0.2);
  }
  boundary ({phi});
  fractions (phi, cs, fs);

I was expecting a circular solid wall but it doesn't seem to be here. What did I do wrong?

Thank you in advance for your help,

Ed

[Edo Frederix]

Dear Pepe,

Regarding your work on embed/metric/two-phase, first, thanks a lot for this effort. It has helped a lot.

In one of your messages you talk about viscosity-embed.h. Have you been able to develop a viscosity.h which is able to deal with embedded boundaries?

The viscosity-embed.h seems to approximate the strain rate tensor as 1/2*dui/dj, instead of the full 1/2*(dui/dj + duj/di). This approximation is valid in the NS equations, in case of incompressibility and constant viscosity. However, in two-phase the viscosity is not constant. Thus, viscosity-embed.h in combination with two-phase will make an error in the viscous term -- not sure how significant it will be. The term that will be missing in the NS equations is dmu/dj*duj/di.

Could someone explain why viscosity-embed.h makes this approximation?

Cheers,

Edo

[jmlopez...@gmail.com]

Hi Edo,

El miércoles, 29 de septiembre de 2021 a la(s) 14:07:59 UTC+2, edofr...@gmail.com escribió:

Dear Pepe,

Regarding your work on embed/metric/two-phase, first, thanks a lot for this effort. It has helped a lot.

Glad it serves to you!

In one of your messages you talk about viscosity-embed.h. Have you been able to develop a viscosity.h which is able to deal with embedded boundaries?

I did a try which was further much improved by Arthur (Ghigo). All our trials are available in te sandbox

Take a look to 

http://basilisk.fr/sandbox/ghigo/src/myviscosity.h

http://basilisk.fr/sandbox/ghigo/src/myembed.h

and as test you could work on

http://basilisk.fr/sandbox/ghigo/src/test-stokes/couette.c

Although our approach works there is still issues with the convergence which is not second order as it should be. So there is stilll works to do.

The viscosity-embed.h seems to approximate the strain rate tensor as 1/2*dui/dj, instead of the full 1/2*(dui/dj + duj/di). This approximation is valid in the NS equations, in case of incompressibility and constant viscosity. However, in two-phase the viscosity is not constant. Thus, viscosity-embed.h in combination with two-phase will make an error in the viscous term -- not sure how significant it will be. The term that will be missing in the NS equations is dmu/dj*duj/di.

 

The only term dismissed (not considered) in viscosity_embed.h is

\[

\nabla \mu \cdot \nabla \mathbf{v} = \partial \mu/ \partial x_i \partial v_i / \partial x_j \mathbf{e}_j

\]

So if the viscosity is uniform, viscosity_embed.h is perfect for an incompressible movement.  A (in principle easy) workaround to generailize viscosity_embed.h is add the aforementioned term explicitly.

Cheers

Pepe

[Edo Frederix]

Thanks for pointing me to those files. Yes, putting the missing term explicitly makes sense too. I added this code

    // In case of embedding, the strain rate tensor is approximated as Sij =

    // 1/2*dui/dj. When having variable viscosity, a term 1/rho*dmu/dj*duj/di is

    // then missing from the momentum equation. This term is added here

    // explicitly.

    #if EMBED

    foreach() {

      foreach_dimension() {

        u.x[] +=

          dt/(2.0*cm[]*rho[]*sq(Delta)+SEPS)

        * (

            (mu.x[1] - mu.x[])*(fm.x[1]*u.x[1] - fm.x[]*u.x[-1])

            #if dimension > 1

          + (mu.y[0,1] - mu.y[])*(fm.x[1]*u.y[1] - fm.x[]*u.y[-1])

            #endif

            #if dimension > 2

          + (mu.z[0,0,1] - mu.z[])*(fm.x[1]*u.z[1] - fm.x[]*u.z[-1])

            #endif

          );

      }

    }

    #endif

after

    mgu = viscosity (u, mu, rho, dt, mgu.nrelax);

Maybe it's helpful to others. Seems to work ok but didn't test it thoroughly. Not sure about stability of this. The metrics are rather approximate.

Cheers,

Edo

[Francesco Picella]

Dear all,

Here you'll find a working example of variable viscosity + embed:

• http://basilisk.fr/sandbox/fpicella/test_embed_variable_viscosity/free_cylinder_in_plane_poiseuille.c
• http://basilisk.fr/sandbox/fpicella/test_embed_variable_viscosity/free_cylinder_in_plane_couette.c

I've essentially updated Arthur's old myviscosity.h

• http://basilisk.fr/sandbox/ghigo/src/myviscosity.h

so that is compatible with the present, stable version of Basilisk:

• http://basilisk.fr/sandbox/fpicella/src/viscosity-embed_variable_viscosity.h

It does it using a fully implicit approach;

still, computation of each velocity component is fully decoupled.

A possible improvement could be working a 2x2 (2D) or 3x3 (3D) system for each cell,

but this goes beyond my scope.

I've tested it in Stokes regime, but I do not see why it shouldn't work at higher Reynolds as well.

I've not tested it using vof...but I let others do it ;)

Hope this code will be useful to someone else!

Best,

Francesco

[Andrés Castillo-Castellanos]

Hi, 

I am interested in using vof + embed. I like Francesco's (and others) implementation, but I have a doubt about the discretization used for the derivatives. 

When we evaluate the usual derivatives, for instance dudx at the x-face, in the cut-cells, we weight the derivatives using fs, like in the diagram below. 

However, when evaluating the "mixed" derivates, for instance the same dudx at the y-face, i feel that we should do something similar. The question is how? The way I see it, we are using u[-1,1] and u[1,0] to compute a gradient (top red arrow, at cell center), and we use u[1,-1] and u[-1,0] to compute another gradient (bottom red arrow at cell center). Then, we average them, to obtain the gradient at the cell face (middle red arrow). For now, both have equal weight, but maybe there is another to weight both gradients? Any suggestions or ideas or a reference that might help?  

Cheers,