<div dir="ltr"><div><div><div><div><div><div>Hi Sang,<br><br></div>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. <br>
<br></div>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.<br>
<br></div>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.<br>
<br></div>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....<br>
<br><br></div>Cheers,<br></div> Dave<br><div><div><div><br><div><div><div><br><br></div></div></div></div></div></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On 17 April 2014 03:11, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div class="">On Wed, Apr 16, 2014 at 5:38 PM, Sang pham van <span dir="ltr"><<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Yes, It does converge toward 0-flow when I refine the mesh !<div><br></div><div>In the code, coordinates are set by </div>
<div> 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);<br>
</div><div>Does it set coordinate for elements center here?</div><div>When I want to try non-uniform grid, should I just modify coordinates attached with the DM, or I need more implementation?</div></div></blockquote><div>
<br></div></div><div>Yes, I think that should be enough.</div><div><br></div><div> Matt</div><div><div class="h5"><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div dir="ltr"><div>
Many thanks.</div>
<div>Sang</div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Apr 16, 2014 at 6:23 PM, Matthew Knepley <span dir="ltr"><<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div>On Wed, Apr 16, 2014 at 4:48 PM, Sang pham van <span dir="ltr"><<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Hi Jed,<div><br></div><div>I modified the ex43 code to enforce no-slip BCs on all boundaries. </div><div>
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). </div>
<div>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?</div></div></blockquote>
<div><br></div></div><div>It sounds like it is due to discretization error. Your incompressibility constraint is not verified element-wise</div><div>(I think ex43 is penalized Q1-Q1), so you can have some flow here. Refine it and see if it converges toward</div>
<div>0 flow.</div><div><br></div><div> Matt</div><div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">
<div>Thank you.</div><div><br></div><div>Sang</div><div><br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Mar 26, 2014 at 2:38 PM, Jed Brown <span dir="ltr"><<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div>Sang pham van <<a href="mailto:pvsang002@gmail.com" target="_blank">pvsang002@gmail.com</a>> writes:<br>
<br>
> Hi Dave,<br>
> I guess you are the one contributed the ex42 in KSP's examples. I want to<br>
> modify the example to solve for stokes flow driven by volume force in 3D<br>
> duct. Please help me to understand the code by answering the following<br>
> questions:<br>
><br>
> 1. Firstly, just for confirmation, the equations you're solving are:<br>
> \nu * \nabla \cdot \nabla U - \nabla P = 0 and<br>
<br>
</div>For variable viscosity, it must be formulated as in the example:<br>
<br>
\nabla\cdot (\nu D U) - \nabla P = 0<br>
<br>
where D U = (\nabla U + (\nabla U)^T)/2<br>
<div><br>
> \nabla \cdot U = 0<br>
><br>
> where U = (Ux,Uy,Uz), \nu is variable viscosity?<br>
><br>
> 2. Are U and P defined at all nodes? (I googled the Q1Q1 element, it looks<br>
> like a box element with U and P defined at 8 corners).<br>
<br>
</div>Yes.<br>
<div><br>
> 3. Are nodes' coordinate defined though the DA coordinates?<br>
<br>
</div>Yes, though they are set to be uniform.<br>
<div><br>
> 4. How can I enforce noslip BC, and where should I plug in volume force?<br>
<br>
</div>Enforce the Dirichlet condition for the entire node.<br>
</blockquote></div><br></div>
</blockquote></div></div><span><font color="#888888"><br><br clear="all"><span><font color="#888888"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</font></span></font></span></div></div>
</blockquote></div><br></div>
</blockquote></div></div></div><div><div class="h5"><br><br clear="all"><div><br></div>-- <br>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
-- Norbert Wiener
</div></div></div></div>
</blockquote></div><br></div>