[petsc-users] distribute and cells mapping.
Olivier Bonnefon
olivier.bonnefon at avignon.inra.fr
Thu Aug 29 04:28:58 CDT 2013
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
> <mailto: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
>>> <mailto: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)
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
>>> <tel:%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 <tel:%2B33%20%280%294%2032%2072%2021%2058>
>
>
> --
> Olivier Bonnefon
> INRA PACA-Avignon, Unité BioSP
> Tel:+33 (0)4 32 72 21 58 <tel:%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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130829/0e0a107d/attachment-0001.html>
More information about the petsc-users
mailing list