[petsc-users] Reading data from a file into a DM created vector

[Matthew Knepley]

On Sat, Nov 1, 2014 at 12:20 PM, Justin Chang <jychang48 at gmail.com> wrote:

> Thanks for the response. On a related note,
>
> In the source code of SNES ex12, where is the PetscSection actually being
> created? I can't seem to find anywhere in the code or its routines where
> the creation has been explicitly called. Because when the
> discretization/problem is define, I imagine a PetscSection for the fields
> and constraints has to be invoked somewhere. It seems to me as soon as
> DMPlexAddBoundary() is called it creates the PetscSectionConstraints but I
> can't seem to really confirm this within the source of that function. The
> reason I want to know this is because if I really were to declare my own
> constraints, I want to make sure that the PetscSection for the field
> variables has already been set up.
>

Much of the mechanism has been moved into the library. The reason I did
this was to make things like
nonlinear multigrid automatic, rather than making the user define a
different section on each level. Let
me explain what is happening. When we call DMGetDefaultSection(), if the
section is not present it
will call the createdefaultsection() member function:


https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/interface/dm.c?at=master#cl-2977

which for Plex is here


https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-5601

It uses the discretization information from PetscDS and the BC information
from DMPlexAddBoundary() to build the inputs for
DMPlexCreateSection()


https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-5739

I realize that call is still limited. For example, it constrains all dofs
on a point, but sometimes you only want to constrain some. In
the most general case, you would just build the Section manually, using the
same tools I use in DMPlexCreateSection():


https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plex.c?at=master#cl-3204

  Thanks,

     Matt


> Thanks,
> Justin
>
>
> Sent from my iPhone
>
> On Nov 1, 2014, at 11:40 AM, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Fri, Oct 31, 2014 at 4:15 PM, Justin Chang <jychang48 at gmail.com> wrote:
>
>> Matt,
>>
>> Thanks for the response. One more question:
>>
>> If I go with approach 2, will manually setting the constraints/indices
>> and its values be compatible with the DMPlexSNESComputeResidual/JacobianFEM
>> routines? When I look at the source code of those routines it seems the
>> constrained values are added to the local solution vectors via
>> DMPlexInsertBoundaryValuesFEM. If I choose not to use DMPlexAddBoundary and
>> wish to manually declare my boundary values, i assume I should call
>> DMPlexProjectField with mode INSERT_BC_VALUES?
>>
>
> Yes, exactly. I use DMPlexProjectFunctionLabelLocal(),
>
>
> https://bitbucket.org/petsc/petsc/src/ea1a9a653238fa98fb68e87b49145608ab6e5301/src/dm/impls/plex/plexfem.c?at=master#cl-576
>
>   Thanks,
>
>      Matt
>
>
>> Thanks,
>> Justin
>>
>> Sent from my iPhone
>>
>> On Oct 30, 2014, at 4:24 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>
>> On Thu, Oct 30, 2014 at 4:16 PM, Justin Chang <jychang48 at gmail.com>
>> wrote:
>>
>>> Matt, thanks for the quick response.
>>>
>>> What about for (dirichlet) boundary conditions? Would it be possible to
>>> do something similar for those, like using those PetscSectionSetConstraint
>>> functions?
>>>
>>
>> Yes. There are generally two ways of handling Dirichlet conditions:
>>
>>   1) Replace those rows of the Jacobian with the identity and put the
>> boundary value in the rhs. I find
>>       this cumbersome for nonlinear solves.
>>
>>   2) Remove these unknowns from the system using
>> PetscSectionSetConstraintDof() and ConstratinIndices().
>>
>> The DMPlexAddBoundary() functions are automating this processing by
>> marking boundaries using DMLabels,
>> and then constraining dofs on those boundaries.
>>
>>   Thanks,
>>
>>     Matt
>>
>>
>>> Thanks,
>>> Justin
>>>
>>> On Thu, Oct 30, 2014 at 3:35 PM, Matthew Knepley <knepley at gmail.com>
>>> wrote:
>>>
>>>> On Thu, Oct 30, 2014 at 3:29 PM, Justin Chang <jychang48 at gmail.com>
>>>> wrote:
>>>>
>>>>> Hi all,
>>>>>
>>>>> So I am writing an FEM code where it reads input data (e.g., auxiliary
>>>>> coefficients, source terms, etc) from a text file. I have preprocessed the
>>>>> data text file so that each vertex point has its corresponding data. For
>>>>> instance, if my source term for a diffusion problem has a sin or cos
>>>>> function of the coordinates, then this data is already tabulated and simply
>>>>> needs to be fed into my PETSc program. The data text file containing both
>>>>> the auxiliary data and the coordinate/connectivities will also be used to
>>>>> provide the input arguments for the DMPlexCreateFromDAG() function.
>>>>>
>>>>> Normally, I would hard-code the sin/cos functions into the source code
>>>>> itself, but i want my program to be "universal" in that it can take as
>>>>> input any user-defined function for not only these auxiliaries but for
>>>>> other things like boundary conditions. I see that PETSc has a lot of
>>>>> examples on how to read data into vectors and matrices, but I guess my
>>>>> question is is there a way to project data from a text file into a vector
>>>>> that has already been created with a defined DM structure?
>>>>>
>>>>
>>>> If you have the functions available, I think it is far more universal
>>>> to use the function itself, since then you can be independent
>>>> of mesh and discretization when specifying input and BC.
>>>>
>>>> However, if you want to read it in, and you guarantee that it matches
>>>> the mesh and discretization, I think the easiest thing to do is
>>>> demand that it come in the same order as the vertices and use
>>>>
>>>> VecGetArray(V, &a);
>>>> for (v = vStart, i = 0; v < vEnd; ++v) {
>>>>   PetscSectionGetDof(s, v, &dof);
>>>>   PetscSectionGetOffset(s, v, &off);
>>>>   for (d = 0; d < dof; ++d) a[off+d] = text_data[i++];
>>>> }
>>>> VecRestoreArray(V, &a);
>>>>
>>>>    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
>
>


-- 
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/20141101/66d87a89/attachment-0001.html>