[petsc-users] DMPlex and Tao

[Matthew Knepley]

On Fri, Aug 1, 2014 at 2:38 PM, Justin Chang <jychang48 at gmail.com> wrote:

> As far as I know this equality constraint should exist outside of the
> large square matrix.
>
> What I was thinking of doing is manually creating an *Nele* by *Ndof*
> matrix, each row will have the same number of non-zeros. Since the size of
> *Ndof* is simply the global number of free dofs formulated by the
> PetscSection, I want to know what's the best way to map the free dofs into
> the global column indices [0 *Ndof*)
>
> For example, in src/dm/impls/plex/examples/tutorials/ex1.c where a box
> mesh with three fields is created, there are a total of 41 dofs to which 8
> of them are constrained. If no constraints were set, then the size of u
> (i.e. *Ndof*) would be 41, but with constraints set it is 33. When you
> view the constrained global vector u, it lists all the dofs in the order as
> specified by the PetscSection but skips all the constrained ones. Would I
> have to use something like DMSet/GetDefaultGlobalSection to do the mapping?
>

I do not understand what you are asking. However, you can just use the
PetscSection for the Ndof field to give you
the local sizes for your matrix.

   Matt


> Thanks,
> Justin
>
>
> On Fri, Aug 1, 2014 at 7:26 AM, Matthew Knepley <knepley at gmail.com> wrote:
>
>> On Wed, Jul 30, 2014 at 5:01 PM, Justin Chang <jychang48 at gmail.com>
>> wrote:
>>
>>> That's good to know thanks. I have another question somewhat related.
>>>
>>> In that quadratic programming problem I have shown earlier, my K matrix
>>> is of global size *Ndof* by *Ndof* where *Ndof* is number of free dofs
>>> (for instance, if I have a 2x2 quadrilateral mesh and want to evaluate
>>> concentration at all the verticies, I will have nine total dofs but if I
>>> have Dirichlet constraints all along the boundary, only my center vertex is
>>> free so *Ndof* = 1. Thus when I simply call DMCreateMatrix(...) it
>>> would give me a 1x1 matrix.)
>>>
>>> However, I want my M constraint matrix to be of size *Nele* by *Ndof*
>>> where *Nele* is the total number of elements. How would I create this
>>> matrix? Because ultimately, I am subjecting all of my free u's to *Nele*
>>> number of equality constraints, and if possible I want this M matrix to be
>>> compatible with the layout from Vec u.
>>>
>>
>> This is not hard in principle. I used to have code that did this, but I
>> took it out during the rewrite because no one
>> was using it and it was complicated. It just generalizes the code in
>> DMPlexPreallocateOperator() to take two
>> PetscSections instead of one.
>>
>> Are you sure that this does not belong inside a large square matrix for
>> the entire problem?
>>
>>   Thanks,
>>
>>      Matt
>>
>>
>>>
>>> On Wed, Jul 30, 2014 at 2:03 PM, Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Wed, Jul 30, 2014 at 2:53 PM, Justin Chang <jychang48 at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi,
>>>>>
>>>>> This might be a silly question, but I am developing a finite element
>>>>> code to solve the non-negative diffusion equation. I am using the DMPlex
>>>>> structure and plan on subjecting the problem to a bounded constraint
>>>>> optimization to ensure non-negative concentrations. Say for instance, I
>>>>> want to solve:
>>>>>
>>>>> min || K*u - F ||^2
>>>>> subject to u >= 0
>>>>>              M*u = 0
>>>>>
>>>>> Where K is the Jacobian, F the residual, u the concentration, and M
>>>>> some equality constraint matrix. I know how to do this via Tao, but my
>>>>> question is can Tao handle Mats and Vecs created via DMPlex? Because from
>>>>> the SNES and TS examples I have seen in ex12, ex62, and ex11, it seems
>>>>> there are solvers tailored to handle DMPlex created data structures.
>>>>>
>>>>
>>>> Yes, TAO can handle objects created by DMPlex since they are just Vec
>>>> and Mat structures. The additions
>>>> to ex12, ex62, and ex11 concern the callbacks from the solver to the
>>>> constructions routines for residual and
>>>> Jacobian. Right now, none of the callbacks are automatic for TAO so you
>>>> would have to code them yourself,
>>>> probably by just calling the routines from SNES or TS so it should not
>>>> be that hard. We are working to get
>>>> everything integrated as it is in the other solvers.
>>>>
>>>>   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
>>
>
>


-- 
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/20140801/5ae207c6/attachment.html>