[petsc-users] ex42 question

[Dave May]

Hi Sang,

The only assumptions made in ex42 are about the mesh topoology, e.g. it
assumes the mesh can be defined via an IJK type representation. Assembly of
the weak form makes no assumptions about elements being affine, or elements
having uniform spacing. The quadrature used is also sufficiently accurate
that you can deform the mesh and the weak form will be correctly
integrated.

I would caution you about using this type of element Q1Q1 for any science
applications which are driven purely by buoyancy variations (e.g. dirichlet
and neumann boundary conditions are of zero value). This element uses
Bochev stabilization and exhibits an error in div(u) which is O(h). The
constant associated with this error can be huge. Also, the constant is very
much dependent on the type of boundary conditions you use.

Actually, one of the tests which shows how careful you should be is similar
to what you tried. Apply u_i n_i=0 on all laterial walls and the base.
Apply tau_{ij} t_j = 0 (no normal flow, zero shear stress on the walls) on
the lateral sides and base. Set \sigma_{ij} =  0 on the surface. e.g. you
want to model flow in a cup of water - the solution should be exactly zero.
As Matt pointed out, this element doesn't have element wise conservative
properties. Run this experiment, and you will see the fluid "compacts" at
the base the domain. Increase the density of the fluid and you will see the
"compaction" increases in size. Add a jump horizontal layering in the
viscosity and you will see "compaction" at the interface.

All these compaction like flow features are local and completely associated
with poor enforcement of div(u). If you refine the mesh, the div(u) error
decreases and the flow field approaches something incompressible. I
strongly encourage you to always monitor div(u) in every element and
consider whether this value is appropriate given your science application.
If its not appropriate, you need to refine the mesh. I gave up with this
element as to obtain meaningful solution I determined I would need 1500^3
elements....


Cheers,
  Dave





On 17 April 2014 03:11, Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, Apr 16, 2014 at 5:38 PM, Sang pham van <pvsang002 at gmail.com>wrote:
>
>> Yes, It does converge toward 0-flow when I refine the mesh !
>>
>> In the code, coordinates are set by
>>   ierr =
>> DMDASetUniformCoordinates(da_prop,0.0+0.5*dx,1.0-0.5*dx,0.0+0.5*dy,1.0-0.5*dy,0.,0);CHKERRQ(ierr);
>> Does it set coordinate for elements center here?
>> When I want to try non-uniform grid, should I just modify coordinates
>> attached with the DM, or I need more implementation?
>>
>
> Yes, I think that should be enough.
>
>    Matt
>
>
>> Many thanks.
>> Sang
>>
>>
>>
>> On Wed, Apr 16, 2014 at 6:23 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>
>>> On Wed, Apr 16, 2014 at 4:48 PM, Sang pham van <pvsang002 at gmail.com>wrote:
>>>
>>>> Hi Jed,
>>>>
>>>> I modified the ex43 code to enforce no-slip BCs on all boundaries.
>>>> I run the code with volume force (0,-1) and isoviscosity. The expected
>>>> result is Vx = Vy = 0 everywhere, and linearly decreasing pressure (from to
>>>> to bottom).
>>>> In the attached is plot of velocity field and pressure, so there is
>>>> still a (light) flow in middle of the domain. Do you know why the solution
>>>> is that, and what should I do to get the expected result?
>>>>
>>>
>>> It sounds like it is due to discretization error. Your incompressibility
>>> constraint is not verified element-wise
>>> (I think ex43 is penalized Q1-Q1), so you can have some flow here.
>>> Refine it and see if it converges toward
>>> 0 flow.
>>>
>>>    Matt
>>>
>>>
>>>> Thank you.
>>>>
>>>> Sang
>>>>
>>>>
>>>>
>>>>
>>>> On Wed, Mar 26, 2014 at 2:38 PM, Jed Brown <jed at jedbrown.org> wrote:
>>>>
>>>>> Sang pham van <pvsang002 at gmail.com> writes:
>>>>>
>>>>> > Hi Dave,
>>>>> > I guess you are the one contributed the ex42 in KSP's examples. I
>>>>> want to
>>>>> > modify the example to solve for stokes flow driven by volume force
>>>>> in 3D
>>>>> > duct. Please help me to understand the code by answering the
>>>>> following
>>>>> > questions:
>>>>> >
>>>>> > 1. Firstly, just for confirmation, the equations you're solving are:
>>>>> >      \nu * \nabla \cdot \nabla U - \nabla P = 0 and
>>>>>
>>>>> For variable viscosity, it must be formulated as in the example:
>>>>>
>>>>>   \nabla\cdot (\nu D U) - \nabla P = 0
>>>>>
>>>>>   where  D U = (\nabla U + (\nabla U)^T)/2
>>>>>
>>>>> >      \nabla \cdot U = 0
>>>>> >
>>>>> >       where U = (Ux,Uy,Uz), \nu is variable viscosity?
>>>>> >
>>>>> >  2. Are U and P defined at all nodes? (I googled the Q1Q1 element,
>>>>> it looks
>>>>> > like a box element with U and P defined at 8 corners).
>>>>>
>>>>> Yes.
>>>>>
>>>>> >  3. Are nodes' coordinate defined though the DA coordinates?
>>>>>
>>>>> Yes, though they are set to be uniform.
>>>>>
>>>>> >  4. How can I enforce noslip BC, and where should I plug in volume
>>>>> force?
>>>>>
>>>>> Enforce the Dirichlet condition for the entire node.
>>>>>
>>>>
>>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20140417/fb22f35d/attachment.html>