[petsc-users] Using DMPlex for Darcys law

[Matthew Knepley]

On Wed, Aug 6, 2014 at 10:47 PM, Justin Chang <jychang48 at gmail.com> wrote:

> Thank you for the answer, i have a few more questions:
>
> 1) In ex62 regarding the functions f0/f1, what's the difference between
> them? Do f0 correspond with test functions and f1 with the derivatives of
> the test functions?
>

The residual F is calculated using

    F(u) = \int_\Omega phi f_0(u, grad u) + grad phi . f_1(u, grad u)

This formulation is in my slides and also this paper
http://arxiv.org/abs/1309.1204.


> 2) For darcy's law, the velocity boundary conditions for each surface is
> simply the normal. i.e. u dot n over Gamma_u = <some value>. In
> DMPlexAddBoundary it seems to me that it constraints every dof
> corresponding to the mesh point (and field), but is there a way where given
> a field (say in this case velocity) I only want to prescribe the x velocity
> on the left and right faces and y velocity on the top and bottom?
>

You are describing a Neumann condition. If you look at ex12, I show both
Dirichlet and Neumann conditions. The Neumann conditions end up
integrating your <some value> function over the boundary (that is what
comes out of integration by parts). So you do that AddBoundary() setting
isEssential to PETSC_FALSE.


> 3) Is there a way/function for DMPlex to output multiple fields to VTK
> files? Or would the field solutions have to be manually separated out from
> the original global solution?
>

VTK is very limited. I would switch to using HDF5. The fields are
automatically split up and there is a Python
script for generating Xdmf files which make the HDF5 usable from Paraview
just like VTK.

If you do not like that, you can try Jed's VTU but I have not tried it yet.

  Thanks,

     Matt


> Thanks,
> Justin
>
>
> On Wed, Aug 6, 2014 at 6:27 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Wed, Aug 6, 2014 at 5:55 PM, Justin Chang <jychang48 at gmail.com> wrote:
>>
>>> Hi,
>>>
>>> So in ex62 of SNES it solves the stokes equation with P2-P1 elements. I
>>> want to solve the darcy equation with P2-P1 elements in a similar fashion.
>>>
>>> The difference between stokes and darcy is that in the weak form of
>>> darcy, the velocity-velocity term is v*u whereas in stokes it's
>>> grad[v]*grad[u]
>>>
>>> In ex62, the function g3_uu (lines 172-183) formulates the
>>> grad[v]*grad[u] term but how would I do it for Darcy's equation where its
>>> simply v*u (assuming drag coefficient is scalar quantity of 1.0)?
>>>
>>
>> So you would want to define g0 instead of g3, since its the u-u term:
>>
>> void g0_uu(const PetscScalar u[], const PetscScalar u_t[], const
>> PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const
>> PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
>> {
>>   const PetscInt Ncomp = spatialDim;
>>   PetscInt       compI;
>>
>>   for (compI = 0; compI < Ncomp; ++compI) g0[compI] = u[compI];
>> }
>>
>> and then register it in SetupProblem:
>>
>>   ierr = PetscDSSetJacobian(prob, 0, 0, g0_uu, NULL,  NULL,
>>  NULL);CHKERRQ(ierr);
>>
>> and do not forget to change the residual definition in f0/f1.
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>>  Thanks,
>>> Justin
>>>
>>
>>
>>
>> --
>> 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/20140808/7f8009a6/attachment.html>