[petsc-users] distribute and cells mapping.
Matthew Knepley
knepley at gmail.com
Fri Aug 30 08:47:39 CDT 2013
On Thu, Aug 29, 2013 at 4:28 AM, Olivier Bonnefon <
olivier.bonnefon at avignon.inra.fr> wrote:
> On 08/28/2013 06:08 PM, Matthew Knepley wrote:
>
> On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon <
> olivier.bonnefon at avignon.inra.fr> wrote:
>
>> On 08/28/2013 10:28 AM, Olivier Bonnefon wrote:
>>
>> On 08/23/2013 04:42 PM, Matthew Knepley wrote:
>>
>> On Fri, Aug 23, 2013 at 9:35 AM, Olivier Bonnefon <
>> olivier.bonnefon at avignon.inra.fr> wrote:
>>
>>> Hello,
>>>
>>> Thanks for your answers, I'm now able to import and distribute a mesh:
>>>
>>
>> You might simplify this to
>>
>> if (rank) {obNbCells = 0; obNbVertex = 0;}
>> ierr =
>> DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>>
>>
>>> if (!rank){
>>> ierr =
>>> DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>>> for (i=0;i<obNbBound;i++){
>>> ierr =DMPlexSetLabelValue(*dm, "marker",
>>> obBoundary[i]+obNbCells, 1);CHKERRQ(ierr);
>>> }
>>> }else {
>>> ierr =
>>> DMPlexCreateFromCellList(comm,dim,0,0,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>>> }
>>>
>>> ierr = DMPlexDistribute(*dm, partitioner, 0,
>>> &distributedMesh);CHKERRQ(ierr);
>>> if (distributedMesh) {
>>> ierr = DMDestroy(dm);CHKERRQ(ierr);
>>> *dm = distributedMesh;
>>> }
>>>
>>> Is it possible to known the resulting partition ? ie, What is the
>>> mapping between the initial cell number and the local cell (used in
>>> DMPlexComputeResidualFEM)?
>>> I need this to write an efficient implementation of the FEM struct
>>> functions f0 and g0, space depending.
>>>
>>
>> Yes, but I really do not think you want to do things that way. I am
>> assuming you want different material models or something
>> in different places. The way I envision that is using a DMLabel to mark
>> up parts of the domain. All labels are automatically
>> distributed with the mesh. Is that what you want?
>>
>> Hello,
>>
>> It is exactly what I need: I'm mobilized a landscape, and the parameters
>> of the model depend of the type of crop. Therefore, I have created a label
>> for each type of crop and I have labeled each triangle with the
>> corresponding label:
>>
>> for (i=0;i<obNbCells;i++){
>> if (labelCells[i]==1){
>> ierr =DMPlexSetLabelValue(*dm, "marker1", i, 1);CHKERRQ(ierr);
>> }else{
>> ierr =DMPlexSetLabelValue(*dm, "marker2", i, 1);CHKERRQ(ierr);
>> }
>> }
>>
>> So, I'm able to mark the triangles, but I'm not able to get this label in
>> the plugin "fem.f0Funcs" and "fem.g0Funcs": These plugins are called by
>> looping on the triangles in the function "FEMIntegrateResidualBatch", but
>> the dm is not available, so I can't use the functions DMPlexGetLabel,
>> DMLabelGetStratumSize and DMLabelGetStratumIS. What is the good way to get
>> the labels in the user plugins of the fem struct ?
>>
>>
> So lets start with the abstract problem so that I can see exactly what
> you want to do. In ex12 (or ex62, etc.) I have a single
> equation, so I do a loop over all cells. This loop takes place in
> DMPlexComputeResidualFEM(). You would instead like
> to do a few loops over sets of cells with different material models, using
> different f0/f1. Is this correct?
>
>
>> Thanks a lot for your help.
>>
>> Olivier B
>>
>> Hello,
>>
>> This is the solution I implemented to get the label level in the plugins
>> "fem.f0Funcs" and "fem.g0Funcs":
>>
>> I need the DM and the index element, so i do:
>> 1) I add some static variables:
>> static DM * spDM[128];
>> static int scurElem[128];
>>
>
> Notice that the DM is available in DMPlexComputeResidualFEM(). Here is
> what the function does:
>
> a) Batches up elements into groups
>
> b) Integrates each group using a call to FEMIntegrateResidualBatch().
> Notice that in 'next' this has
> changed to PetscFEIntegrateResidual() since we have added a few FEM
> classes to make things
> simpler and more flexible.
>
> What you can do, I think, to get what you want is:
>
> a) Write a new MY_DMPlexComputeResidualFEM() to do a few loops. This
> is supplied to your app using
>
> ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode
> (*)(DM,Vec,Vec,void*)) MY_DMPlexComputeResidualFEM, &user);CHKERRQ(ierr);
> ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode
> (*)(DM,Vec,Mat,Mat,MatStructure*,void*)) MY_DMPlexComputeJacobianFEM,
> &user);CHKERRQ(ierr);
>
> just as in the examples. You could use different f0/f1 for each loop
> somehow.
>
> b) Write a new PetscFEIntegrateResidual() that does what you want. The
> easiest way to do this is create
> a new PetscFE subclass, since they only really do one thing which
> is these integrals. I can help you.
>
> HOWEVER, if what you really want to do is get coefficient information into
> f0/f1 instead of a different physical model,
> then you can do something easier that we just put in. You can layout a
> coefficient, like nu in
>
> \div \nu \grad u = \rho
>
> and provide a DM for \nu. This will be passed all the way down inside
> until f0 gets
>
> f0Func(u, gradU, nu, gradNu, x, f0)
>
> so that the pointwise values of the your coefficient and its gradient
> are available to your physics.
>
> I am sure there will be questions about this, but the first thing to do
> is get entirely clear what you want to do.
>
> Hello,
>
> This is the 2D systems I want to simulate:
>
> system 1:
> 0=\div \rho(x) \grad u + r(x)*u(1-u)
>
I have just pushed a version of this system to SNES ex12 (you have to use
the 'next' branch). I see at least
two options:
a) You have rho(x) and r(x) as explicit functions
This is easy since the f0,f1 functions get x explicitly. This is the
-variable_coefficient analytic mode of ex12.
b) You have another FEM field defining rho(x) and r(x)
This is the -variable_coefficient field mode of ex12. You can use an
element different from that for u is defining
these coefficients. For example, I can use P0 for rho and P2 for u, or
P1 for both.
I can try and help you modify that example to solve this problem if you
want, or you can continue to manage the
lower level stuff since that gives you more flexibility.
Thanks,
Matt
> corresponding to the stationary state of the time dependent problem:
> system 2:
> \partial_t u = \div \rho(x) \grad u + r(x)*u(1-u)
>
> It represents the diffusion of a specie throw a 2D space, u is the density
> of this specie.
>
> I want also to study the effect of a predator p:
> system 3:
> \partial_t u = \div \rho_u(x) \grad u + r(x)*u(1-u) - \beta (x) p*u
> \partial_t p = \div \rho_p(x) \grad p - \delta (x) p + \gamma (x) p*u
>
> I'm focused on the (system 1).
>
> About the geometry:
> The geometry come from the landscape composed of crops. There are
> different type of crops. The functions ( \rho(x), r(x), \delta (x), \beta
> (x)) depend on this type of crops.
>
> I'm focused on the (system 1). I'm working from the ex12.c. Therefore, the
> plungins functions are:
>
> f0_u(u,grad_u,x,f0){
> ...
> f0= r(x)*u*(1-u)
> ...
> }
>
> f1_u(u,grad_u,x,f1){
> ...
> f1[comp*dim+d] = rho(x)*gradU[comp*dim+d];
> ...
> }
>
> g0_u(u,grad_u,x,g0){
> ...
> g0=-r(x)*(1-2*u)
> ...
> }
>
> g3_uu(u,grad_u,x,g3){
> ...
> g3[((compI*Ncomp+compI)*dim+d)*dim+d] = \rho(x);
> ...
> }
>
> For an efficient implementation of theses plugins, I have to know the type
> of crop. If I well understand your previous mail, you propose to me to
> defined my own struct PetscFEM adding my useful parameters (crop type for
> example). I have to overload the functions DMPlexComputeResidualFEM,
> DMPlexComputeJacobianFEM, FEMIntegrateResidualBatch,
> FEMIntegrateJacobianActionBatch. I agree, thanks a lot.
>
>
> Regards,
>
> Olivier B
>
>
> Thanks,
>
> Matt
>
> 2) I overload the DMPlexComputeJacobianFEM with :
>> PetscErrorCode MY_DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat
>> JacP, MatStructure *str,void *user)
>> {
>>
>> PetscMPIInt rank;
>> PetscErrorCode ierr;
>>
>> PetscFunctionBeginUser;
>> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
>> spDM[rank]=&dm;
>> PetscFunctionReturn(DMPlexComputeJacobianFEM(dm, X, Jac,JacP, str,user));
>>
>> }
>> 3) overload FEMIntegrateResidualBatch adding code:
>> .
>> .
>> for (e = 0; e < Ne; ++e) {
>> scurElem[rank]=e;//added ligne
>> .
>> .
>>
>> So that, I can get the label level using DMPlexHasLabel and
>> DMLabelGetValue
>>
>> I'm sure this solution is awful, and works only in this version, but i
>> didn't find a better way to get the label in the plugins fem struc. Do you
>> know the correct way to do that ??
>>
>> Thanks,
>>
>> Olivier B
>>
>>
>> Thanks,
>>
>> Matt
>>
>>
>>> Regards,
>>>
>>> Olivier B
>>>
>>> --
>>> Olivier Bonnefon
>>> INRA PACA-Avignon, Unité BioSP
>>> Tel: +33 (0)4 32 72 21 58 <%2B33%20%280%294%2032%2072%2021%2058>
>>>
>>>
>>
>>
>> --
>> 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
>>
>>
>>
>> --
>> Olivier Bonnefon
>> INRA PACA-Avignon, Unité BioSP
>> Tel: +33 (0)4 32 72 21 58
>>
>>
>>
>> --
>> Olivier Bonnefon
>> INRA PACA-Avignon, Unité BioSP
>> Tel: +33 (0)4 32 72 21 58
>>
>>
>
>
> --
> 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
>
>
>
> --
> Olivier Bonnefon
> INRA PACA-Avignon, Unité BioSP
> Tel: +33 (0)4 32 72 21 58
>
>
--
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/20130830/db8bba77/attachment.html>
More information about the petsc-users
mailing list