Different ways of setting free-slip BCs have differing results

[jcgraciosa]

This is the issue I mentioned earlier in the meeting and posting here for documentation.

If I use different ways (shown below) of setting free-slip BCs then I get different results.

Method 1:

stokes.add_dirichlet_bc((0.0,), "Bottom", (1,))
stokes.add_dirichlet_bc((0.0,), "Top", (1,))
stokes.add_dirichlet_bc((0.0,), "Left", (0,))
stokes.add_dirichlet_bc((0.0,), "Right", (0,))

Method 2:

stokes.add_dirichlet_bc((sympy.oo,0.0), "Bottom")
stokes.add_dirichlet_bc((sympy.oo, 0.0), "Top")
stokes.add_dirichlet_bc((0.0,sympy.oo), "Left")
stokes.add_dirichlet_bc((0.0,sympy.oo), "Right")

Results:

Cartesian convection: Resolution = 16: Method 1:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -29.32204257096935
max. v (component): 23.151495282044586

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.873897927155742
max. v (component): 23.067035407567737

Resolution = 32: Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINE_SEARCH iterations 8
min. v (component): -6.316221446999447
max. v (component): 4.949063245632038

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.820669550518556
max. v (component): 23.160314007131472

Resolution = 64 Method 1:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -29.398087782064863
max. v (component): 23.421398167185224

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -27.964331720326253
max. v (component): 22.900686663482414

SolC benchmark: Resolution = 16: Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -0.023157758823695623
max. v (component): 0.02315529270672989

Resolution = 32: Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
min. v (component): -0.023190411986986556
max. v (component): 0.02319011766071459

Resolution = 64: Method 1:

Nonlinear stokes_ solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
min. v (component): -2.738612787525831
max. v (component): 0.8660254037844388

Method 2:

Nonlinear stokes_ solve converged due to CONVERGED_FNORM_RELATIVE iterations 2
min. v (component): -0.023198938179237818
max. v (component): 0.02319896279134057

Additional note:

• used default petsc settings and same stokes solver tolerance value

[bknight1]

The DIVERGED_LINEAR_SOLVE error in SolC suggests the BC isn't applied correctly.

if you try:

stokes.add_dirichlet_bc((0.0,0.0), "Bottom", (1,))
stokes.add_dirichlet_bc((0.0,0.0), "Top", (1,))
stokes.add_dirichlet_bc((0.0,0.0), "Left", (0,))
stokes.add_dirichlet_bc((0.0,0.0), "Right", (0,))

both method 1 and method 2 should produce the same result ( may have to be something like sympy.Matrix([0.,0.]) instead of (0.0,0.0) )

[lmoresi]

As discussed today, we need to choose a better interface and change all the examples to the new one.

[lmoresi]

We discussed switching from

stokes.add_dirichlet_bc((sympy.oo,0.0), "Bottom")

to this:

stokes.add_dirichlet_bc((None,0.0), "Bottom")

But the problem with this choice is that we can't pass in a sympy.Matrix with that form because

sympy.Matrix([None, 1, None]) is deprecated

1: SymPyDeprecationWarning:

non-Expr objects in a Matrix is deprecated. Matrix represents a mathematical matrix. To represent a container of non-numeric entities, Use a list of lists, TableForm, NumPy array, or some other data structure instead.

See https://docs.sympy.org/latest/explanation/active-deprecations.html#deprecated-non-expr-in-matrix for details.

This has been deprecated since SymPy version 1.9. It will be removed in a future version of SymPy.

@julesghub @bknight1 @jcgraciosa - what do you think ?

[bknight1]

Hmmm that's unfortunate but at least we've found out now.

I like sympy.oo under the hood for cancelling stuff out but not a big fan of it being user-facing. To me it's not intuitive that you'd use sympy.oo to have that component unconstrained. I think something like method 1 above is more intuitive, if the component is not set then it's unconstrained. I don't think it's that big of a deal though as once you've seen an example you'd know how to set free slip/no slip boundaries etc

[lmoresi]

I think it would be fine to use the UWGEO style bc description [None, 0, None] and sanitise it for sympy once we receive it in the solver. I like the idea that we are careful to specify the whole boundary in one go as it is very easy to specify inconsistent bc's otherwise.

[bknight1]

So convert it to a sympy matrix under the hood, and covert None to sympy.oo, or something along those lines? I think that also sounds good!

[lmoresi]

Yes, basically. We currently do some processing to

1) make sure it is a sympy.Matrix type 2) determine which components are sympy.oo to construct a components list.

We do ultimately need a matrix for compilation as that's now assumed in the JIT code. We could leave it as sympy.oo because the compiler puts something useful into the C code when it sees that and it might help for debugging.

The real question is how to provide a pretty-printed view for the users. Perhaps we use a sympy symbol with a display form that is the empty set.

[NengLu]

We may also need to include an extended boundary condition setting to make it possible to add extra node index sets (e.g., some internal nodes) flexibly.

In UWGeo style, it is like: set_velocityBCs(left=[0,None], right=[0,None] top=[None,0], bottom=[None,0],nodeSets=[[fn,0],interface_nodes])

Now, it could be done in uw3 by embedding the mesh with an internal surface which is added as a label by the gmsh. stokes.add_natural_bc((0,0), "Internal",(0,))

Could we make it more direct in _petsc_generic_snessolvers. _setup_discretisation() to add nodes in the boundary condition beside the "label"?

[lmoresi]

Have a look at the Kernel examples which do use an internal boundary.

Note that the PETSc interface for bc's is a bit hampered by its history and natural v essential bcs are a little bit different. At the moment we specify the natural bc's through a row matrix for all components because it is not entirely obvious what the components mean (to the user) on these embedded surfaces.

That's why I want to move to the same approach for other bc's.

[lmoresi]

@jcgraciosa - Issue noted, discussion above, can you finalise the proposed interface changes and we'll decide how to implement them.

[jcgraciosa]

As a user, I think it's more intuitive to do the (None, 0., None) style. Also, since this is how BCs are set in UWGeo, sticking to it would be easier for long-time users (or those who need to use both uw3 and UWGeo) to transition.

[julesghub]

Yep - please go with the UWGeo style @jcgraciosa

[lmoresi]

Internally, the sympy Matrix representation will need to use either sympy.oo or sympy.nan because None is not permitted as an entry in a sympy Matrix. I suggest a conversion to sympy.oo because it prints out better if you ask what boundary conditions have been set.

[lmoresi]

@jcgraciosa - close if completed, please !

[lmoresi]

Sorry - not mine to close.