[petsc-users] DMSNESSet<Function,Jacobian>

Matthew Knepley knepley at gmail.com
Tue Apr 2 00:28:52 CDT 2013


On Tue, Apr 2, 2013 at 4:18 PM, Dharmendar Reddy <dharmareddy84 at gmail.com>wrote:

> Hello,
>           I should be able to generate petscFEM struct with my code, i
> have the required data, I just need to prepare it in the format used by
> DMPLexComputeResidualFEM. I will give it a shot.
>

Cool. If its not general enough, you can always write the loop yourself, or
we can look into expanding it.


> What is the pStart and pEnd in the call to DMLabelCreateIndex referring to
> ? Also,my code complains if i use DMLabel :: labelName
>

It builds an index, so it needs bounds on the possible points you will be
mapping.


> This name does not have a type, and must have an explicit type.
>

Ah, I did not put DMLabel into Fortran. Goes on the list.


> For now i can do without DMLabelCreateIndex as fortran interface to
> DMLabel doesn't seem to exist.
>

Yes, its only an optimization.

  Thanks,

     Matt


> On Mon, Apr 1, 2013 at 11:02 PM, Matthew Knepley <knepley at gmail.com>wrote:
>
>> On Tue, Apr 2, 2013 at 1:47 PM, Dharmendar Reddy <dharmareddy84 at gmail.com
>> > wrote:
>>
>>> Hello,
>>>          I looked at the source code of SNESSetFunction after i sent the
>>> email. From what i understood, SNESSetFunction calls DMSNESSetFunction
>>> internally. I was like, why not just use SNESSetFunction with out passing
>>> DM. But i get it now, not all problems have to use DM. If i am using DM i
>>> should rather use DMSNESSetFunction. A better way is DMSNESSetFunctionLocal
>>> right ? DMSNesSetFunction i called inside DMSNESSetFunctionLocal.
>>>
>>>
>>> I did look at DMPlexComputeResidualFEM but that requires userContext to
>>> have PetscFEM object. Is PetscFEM object accisable from Fortran or is it
>>> genrated using
>>>
>>
>> The object is just a struct, and has only information you must have for
>> FEM, namely
>> a quadrature rule, and tabulation of the basis functions are quadrature
>> points.
>>
>>
>>> bin/pythonscripts/PetscGenerateFEMQuadrature.py
>>>
>>> This script just generates a header with the struct information. You can
>> make this header
>> yourself and obviate the need for the script.
>>
>>
>>> I structured my code such that, user provides the Linear and bilinear
>>> form subroutines:
>>>
>>> subroutine biLinearFrom(u,v,cell,Kelm,eqnCtx)
>>> type(FunctionSpace) :: u,v
>>> type(FEMCell) :: cell
>>> PetscScalar :: Kelm
>>>  class(*),pointer :: eqnCtx
>>> end subroutine
>>>
>>> Function space object provide the routines for evaluation of basis
>>> functions and gradients, given cell coordinates, and cell degrees of
>>> freedom.
>>>
>>
>> Can you just use these to create the PetscFEM struct?
>>
>>
>>>
>>> So i should be able to get the cell geometry information and cell degree
>>> of freedom data via DMPlexVecGet/SetClosure right ?
>>>
>>
>> Yes, although tell me why DMPlexComputeCellGeometry() does not work, and
>> if not you can look at the code
>> to see me getting the coordinate information.
>>
>>
>>> At this point i was able to create a DMPlex object from a mesh created
>>> by gmsh. I used DMPLexCreateLabel to add the physical region labels i
>>> created in GMSh . So no i have a mapping from given a region label, list of
>>> cells in that region.
>>>
>>> If i am doing a physical region wise assemble of the problem, i could
>>> use the DMPlexgetStratumIS for each regionlabel.
>>>
>>
>> Yes.
>>
>>
>>>  since i am doing a cell wise assemble, i would be interested in a
>>> mapping from cell id to  regionId (i can have a local array of labels
>>> indexed by region Id) so that i can pass the material information to
>>> FEMCell.  Right now i have a integer array which holds the map from cellId
>>> to regionId.
>>>
>>> I was wondering if i could add a integer section to DM object ?
>>>
>>
>> You use DMLabelGetValue(), and to trade memory for speed
>> call DMLabelCreateIndex() first if you want.
>>
>>
>>>  I also need to have some other auxiliary scalar data associated with
>>> the DM. You suggested in one of the earlier email  that i created another
>>> DM using DMplexClone ? If i clone a object is it creating another copy of
>>> the mesh in memory or is it just storing pointers ?
>>>
>>
>> Just a pointer.
>>
>>
>>>  I was thinking i could create a section using DMPlexCreateSection and
>>> pass that section to DMPlexVecGetClosure when i have to access auxilary
>>> data.
>>>
>>
>> This is exactly what I was suggesting, but you have the section be the
>> "default section" of the cloned DM. That way
>> it interacts with all PETSc functions nicely and you only have the extra
>> memory for a DM object which is tiny.
>>
>>     Matt
>>
>>  Basically the problem is to solve Poisson equation with current
>>> continuity equation for a semiconductor. I solve the problem by solving
>>> each equation until the solutions of both are consistent. For Poisson
>>> problem, the degree freedom is potential and takes charge density as input.
>>> For continuation equation degree of freedom is charge density and take
>>> potential as input.
>>>
>>> Do you have alternative approaches to what i am thinking ?
>>>
>>> On Mon, Apr 1, 2013 at 8:24 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>
>>>> On Tue, Apr 2, 2013 at 10:16 AM, Dharmendar Reddy <
>>>> dharmareddy84 at gmail.com> wrote:
>>>>
>>>>> Hello,
>>>>>         I am confused about the usage of SNESSetFunction,
>>>>> DMSNESSetFunction, DMSNESSetFunctionLocal and the corresponding Jacobin
>>>>> routines.
>>>>>
>>>>> At present my code is set up as follows.
>>>>>
>>>>> I initiate SNES with getFunction and getJacobain routines required as
>>>>> per SNESSet<Function,Jacobian>. I pass a user context  which has the mesh
>>>>> and field layout information. I assemble the residual function and Jacobin
>>>>> element wise using the information from user context. My code at this point
>>>>> is serial, so i run the element loop from 1 to numberofTotalElements.
>>>>>
>>>>> Now i want to switch to using DMPlex. How should i change the
>>>>> interface to my code.
>>>>>
>>>>
>>>> You can always look at SNES ex62
>>>> http://www.mcs.anl.gov/petsc/petsc-dev/src/snes/examples/tutorials/ex62.c.html
>>>>
>>>>
>>>>> If i set the dm for snes using SNESSetDM, i can see that the DM can be
>>>>> accessed via SNESGetDM inside getFunction and getJacobian, then i get the
>>>>> elementStartID and elementEndID to run the loop for assmbler. But now i
>>>>> will have to assemble to local vector into global vector right ?
>>>>>
>>>>
>>>> No, here is the sequence:
>>>>
>>>>   1) SNESSetDM(snes, dm)
>>>>
>>>>   2) DMSNESSetFunctionLocal(dm, userResidual, userCtx)
>>>>
>>>> where we have
>>>>
>>>>   userResidual(DM dm, Vec X, Vec F, void *userCtx)
>>>>
>>>> Both X and F are local vectors which I normally interact with using
>>>> DMPlexVecGet/SetClosure(). If you
>>>> are using FEM, you can DMPlexComputeResidualFEM() here and
>>>> use DMPlexSetFEMIntegration() to
>>>> input point-wise physics functions as is done in ex62.
>>>>
>>>>
>>>>> Looks like the DMSNESSetFunctionLocal will assemble the local
>>>>> vector into global vector. The function provided to this should just
>>>>> evaluate the local vector. Did i understand this right ?
>>>>>
>>>>> I am confused about passing the DM in DMSNESSet<Function,Jacobian>.
>>>>> When the DM can be accessed via SNESGetDM, why do we pass it again
>>>>> explicitly ?
>>>>>
>>>>
>>>> You are not passing the SNES, so where would it come from in this call?
>>>>
>>>
>>>
>>>
>>>
>>>
>>>>
>>>>    Matt
>>>>
>>>>
>>>>>
>>>>> --
>>>>> -----------------------------------------------------
>>>>> Dharmendar Reddy Palle
>>>>> Graduate Student
>>>>> Microelectronics Research center,
>>>>> University of Texas at Austin,
>>>>> 10100 Burnet Road, Bldg. 160
>>>>> MER 2.608F, TX 78758-4445
>>>>> e-mail: dharmareddy84 at gmail.com
>>>>> Phone: +1-512-350-9082
>>>>> United States of America.
>>>>> Homepage: https://webspace.utexas.edu/~dpr342
>>>>>
>>>>
>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>
>>>
>>>
>>> --
>>> -----------------------------------------------------
>>> Dharmendar Reddy Palle
>>> Graduate Student
>>> Microelectronics Research center,
>>> University of Texas at Austin,
>>> 10100 Burnet Road, Bldg. 160
>>> MER 2.608F, TX 78758-4445
>>> e-mail: dharmareddy84 at gmail.com
>>> Phone: +1-512-350-9082
>>> United States of America.
>>> Homepage: https://webspace.utexas.edu/~dpr342
>>>
>>
>>
>>
>> --
>> 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
>>
>
>
>
> --
> -----------------------------------------------------
> Dharmendar Reddy Palle
> Graduate Student
> Microelectronics Research center,
> University of Texas at Austin,
> 10100 Burnet Road, Bldg. 160
> MER 2.608F, TX 78758-4445
> e-mail: dharmareddy84 at gmail.com
> Phone: +1-512-350-9082
> United States of America.
> Homepage: https://webspace.utexas.edu/~dpr342
>



-- 
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/20130402/a2ed2960/attachment.html>


More information about the petsc-users mailing list