<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class=""><div class=""><br class=""></div> Please open a MR from your fork and we'll take a look at it.<div class=""><br class=""></div><div class=""> Thanks</div><div class=""><br class=""></div><div class=""> Barry</div><div class=""><br class=""><div><br class=""><blockquote type="cite" class=""><div class="">On Feb 21, 2022, at 6:20 PM, Nicolás Barnafi <<a href="mailto:nabw91@gmail.com" class="">nabw91@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><div dir="ltr" class="">Dear all, <div class=""><br class=""></div><div class="">I carried out the required impementation of the PCSetCoordinates for the fieldsplit and wanted to put it upstream. I am not sure about what could be missing (if anything), so perhaps I should open the pull request from a fork in order for you to take a look. Would you advise this or think it is better if I iterate with someone by email? </div><div class=""><br class=""></div><div class="">Best regards, </div><div class="">Nicolas</div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Jan 13, 2022 at 7:21 PM Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> wrote:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class=""><div dir="ltr" class="">On Thu, Jan 13, 2022 at 1:15 PM Nicolás Barnafi <<a href="mailto:nabw91@gmail.com" target="_blank" class="">nabw91@gmail.com</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class="">Dear all, <div class=""><br class=""></div><div class="">I have created a first implementation. For now it must be called after setting the fields, eventually I would like to move it to the setup phase. The implementation seems clean, but it is giving me some memory errors (free() corrupted unsorted chunks).</div><div class=""><br class=""></div><div class="">You may find the code below. After some work with gdb, I found out that the errors appears when calling the ISDestroy(&is_coords) line, which to me is not very clear, as I am indeed within the while scope creating and then destroying the is_coords object. I would greatly appreciate it if you could give me a hint on what the problem is. </div><div class=""><br class=""></div><div class="">After debugging this, and working on your suggestion, I will open a PR. </div><div class=""><br class=""></div><div class="">Best regards, </div><div class="">NB</div><div class=""><br class=""></div><div class=""><br class=""></div><div class="">----- CODE ------</div><div class=""><br class=""></div><div class="">static PetscErrorCode PCSetCoordinates_FieldSplit(PC pc, PetscInt dim, PetscInt nloc, PetscReal coords[])<br class="">{<br class=""> PetscErrorCode ierr;<br class=""> PC_FieldSplit *jac = (PC_FieldSplit*)pc->data;<br class=""> PC_FieldSplitLink ilink_current = jac->head;<br class=""> PC pc_current;<br class=""> PetscInt nmin, nmax, ii, ndofs;<br class=""> PetscInt *owned_dofs; // Indexes owned by this processor<br class=""> PetscReal *coords_block; // Coordinates to be given to the current PC<br class=""> IS is_owned;<br class=""><br class=""> PetscFunctionBegin;<br class=""> // Extract matrix ownership range to then compute subindexes for coordinates. This results in an IS object (is_owned).<br class=""> // TODO: This would be simpler with a general MatGetOwnershipIS (currently supported only by Elemental and BLAS matrices).<br class=""> ierr = MatGetOwnershipRange(pc->mat,&nmin,&nmax);CHKERRQ(ierr);<br class=""> ndofs = nmax - nmin;<br class=""> ierr = PetscMalloc1(ndofs, &owned_dofs); CHKERRQ(ierr);<br class=""> for(PetscInt i=nmin;i<ndofs;++i)<br class=""> owned_dofs[i] = nmin + i;<br class=""> ierr = ISCreateGeneral(MPI_COMM_WORLD, ndofs, owned_dofs, PETSC_OWN_POINTER, &is_owned); CHKERRQ(ierr);<br class=""></div></div></blockquote><div class=""><br class=""></div><div class="">Here you tell PETSc to take control of the memory for owned_dofs, but below you PetscFree(owned_dofs). This is not compatible.</div><div class="">You should just destroy the IS when you are done.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr" class=""><div class=""> // For each IS, embed it to get local coords indces and then set coordinates in the subPC.<br class=""> ii=0;<br class=""> while(ilink_current)<br class=""> {<br class=""> IS is_coords;<br class=""> PetscInt ndofs_block;<br class=""> const PetscInt *block_dofs_enumeration; // Numbering of the dofs relevant to the current block<br class=""><br class=""> ierr = ISEmbed(ilink_current->is, is_owned, PETSC_TRUE, &is_coords); CHKERRQ(ierr); // Setting drop to TRUE, although it should make no difference.<br class=""> ierr = PetscMalloc1(ndofs_block, &coords_block); CHKERRQ(ierr);<br class=""> ierr = ISGetLocalSize(is_coords, &ndofs_block); CHKERRQ(ierr);<br class=""> ierr = ISGetIndices(is_coords, &block_dofs_enumeration); CHKERRQ(ierr);<br class=""><br class=""> // Having the indices computed and the memory allocated, we can copy the relevant coords and set them to the subPC.<br class=""> for(PetscInt dof=0;dof<ndofs_block;++dof)<br class=""> for(PetscInt d=0;d<dim;++d)<br class=""> {<br class=""> coords_block[dim*dof + d] = coords[dim * block_dofs_enumeration[dof] + d];<br class=""> // printf("Dof: %d, Global: %f\n", block_dofs_enumeration[dof], coords[dim * block_dofs_enumeration[dof] + d]);<br class=""> }<br class=""> ierr = ISRestoreIndices(is_coords, &block_dofs_enumeration); CHKERRQ(ierr);<br class=""> ierr = ISDestroy(&is_coords); CHKERRQ(ierr);<br class=""> ierr = KSPGetPC(ilink_current->ksp, &pc_current); CHKERRQ(ierr);<br class=""> ierr = PCSetCoordinates(pc_current, dim, ndofs_block, coords_block); CHKERRQ(ierr);<br class=""> ierr = PetscFree(coords_block); CHKERRQ(ierr);<br class=""> if(!pc_current)<br class=""> SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Setting coordinates to PCFIELDSPLIT but a subPC is null.");<br class=""><br class=""> ilink_current = ilink_current->next;<br class=""> ++ii;<br class=""> }<br class=""> ierr = PetscFree(owned_dofs); CHKERRQ(ierr);<br class=""> PetscFunctionReturn(0);<br class="">}<br class=""></div></div><br class=""><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, Jan 12, 2022 at 6:22 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" class="">bsmith@petsc.dev</a>> wrote:<br class=""></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div class=""><br class=""><div class=""><br class=""><blockquote type="cite" class=""><div class="">On Jan 11, 2022, at 9:51 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank" class="">knepley@gmail.com</a>> wrote:</div><br class=""><div class=""><div dir="ltr" class=""><div dir="ltr" class="">On Tue, Jan 11, 2022 at 3:31 PM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" class="">bsmith@petsc.dev</a>> wrote:<br class=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br class="">
Nicolas,<br class="">
<br class="">
For "simple" PCFIELDSPLIT it is possible to pass down the attached coordinate information. By simple I mean where the splitting is done by fields and not by general lists of IS (where one does not have enough information to know what the coordinates would mean to the subPCS).<br class="">
<br class="">
Look in fieldsplit.c PCFieldSplitSetFields_FieldSplit() where it does the KSPCreate(). I think you can do a KSPGetPC() on that ksp and PCSetCoordinates on that PC to supply the coordinates to the subPC. In the function PCFieldSplitSetIS_FieldSplit() you can also attach the coordinates to the subPCs IF defaultsplit is true.<br class="">
<br class="">
Sadly this is not the full story. The outer PC will not have any coordinates because calling PCSetCoordinates on a PCFIELDSPLIT does nothing since fieldsplit doesn't handle coordinates. So you need to do more, you need to provide a PCSetCoordinates_FieldSplit() that saves the coordinates in new entries in the PC_FieldSplit struct and then in PCFieldSplitSetFields_FieldSplit() you need to access those saved values and pass them into the PCSetCoordinates() that you call on the subPCs. Once you write <br class="">
PCSetCoordinates_FieldSplit() you need to call <br class="">
<br class="">
ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_FieldSplit);CHKERRQ(ierr); <br class="">
<br class="">
inside PCCreate_FieldSplit().<br class="">
<br class="">
Any questions just let us know.<br class=""></blockquote><div class=""><br class=""></div><div class="">I will add "Why is this so cumbersome?". This is a workaround in order to get geometric information into GAMG. It should really be PCGAMGSetCoordinates(), which</div><div class="">are used to calculate the rigid body modes, and assume a bunch of stuff about the coordinate space. This would not help you, because it would still force you to pull</div><div class="">out the correct subPC. The "right" way now to give geometric information to a TS/SNES/KSP/PC is through a DM, which are passed down through PCFIELDSPLIT,</div><div class="">PCPATCH, etc. However they are heavier weight than just some coordinates.</div></div></div></div></blockquote><div class=""><br class=""></div> This is not cumbersome at all. It is a simple natural way to pass around coordinates to PC's and, when possible, their children. </div><div class=""><br class=""></div><div class=""> Barry</div><div class=""><br class=""></div><div class=""> Note that we could also have a DMGetCoordinates() that pulled coordinates from a DM (that happended to have them) in this form associated with the PC and the PC could call it to get the coordinates and use them as needed. But this simple PCSetCoordinates() is a nice complement to that approach.</div><div class=""><br class=""><blockquote type="cite" class=""><div class=""><div dir="ltr" class=""><div class="gmail_quote"><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Barry<br class="">
<br class="">
<br class="">
> On Jan 11, 2022, at 11:58 AM, Nicolás Barnafi <<a href="mailto:nabw91@gmail.com" target="_blank" class="">nabw91@gmail.com</a>> wrote:<br class="">
> <br class="">
> Dear community, <br class="">
> <br class="">
> I am working on a block preconditioner, where one of the blocks uses HYPRE's AMS. As it requires the coordinates of the dofs, I have done so to the PC object. I expected the coordinates to be inherited down to the subblocks, is this not the case? (it seems so as I couldn't find a specialized FIELDSPLIT SetCoordinates function). <br class="">
> <br class="">
> If this feature is missing, please give me some hints on where to add the missing function, I would gladly do it. If not, please let me know why it was dismissed, in order to do things the hard way [as in hard-coded ;)]. <br class="">
> <br class="">
> Kind regards, <br class="">
> Nicolas<br class="">
<br class="">
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></div>
</div></blockquote></div><br class=""></div></blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class="">Nicolás Alejandro Barnafi Wittwer</div></div></div></div>
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class="">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">-- Norbert Wiener</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a><br class=""></div></div></div></div></div></div></div></div>
</blockquote></div><br clear="all" class=""><div class=""><br class=""></div>-- <br class=""><div dir="ltr" class="gmail_signature"><div dir="ltr" class=""><div class=""><div dir="ltr" class="">Nicolás Alejandro Barnafi Wittwer</div></div></div></div>
</div></blockquote></div><br class=""></div></body></html>