[petsc-users] Update DMDA attached to DMSWARM
Matthew Knepley
knepley at gmail.com
Fri Jan 24 07:20:31 CST 2025
On Fri, Jan 24, 2025 at 4:41 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es> wrote:
> Dear Matt, the error was in the implementation of the volume expansion
> function. I updated it, and it works finte under finite domains. However,
> if I include periodic boundary conditions the volume of the cell does not
> accommodate the volume expansion of the particles. The deformation gradient
> is not the identity… I guess I am missing the fine detail on how periodic
> bcc are implemented in DMDA mesh, I’m right?
>
DMDA identifies vertices using a VecScatter to implement periodic BC. This
should be insensitive to coordinates. However, I don't think the algorithm
below is correct for local coordinates. You use GlobalToLocal(), which
means that some global coordinate "wins" for each local cell, so cells on
the periodic boundary can be wrong. I would set the local coordinates by
hand as well.
Thanks,
Matt
> Thanks,
> Miguel
>
> static PetscErrorCode Volumetric_Expansion_DMDA(DM * da,
> const Eigen::Matrix3d& F) {
>
> PetscInt i, j, mstart, m, nstart, n, pstart, p, k;
> Vec local, global;
> DMDACoor3d ***coors, ***coorslocal;
> DM cda;
>
> PetscFunctionBeginUser;
> PetscCall(DMGetCoordinateDM(*da, &cda));
> PetscCall(DMGetCoordinates(*da, &global));
> PetscCall(DMGetCoordinatesLocal(*da, &local));
> PetscCall(DMDAVecGetArray(cda, global, &coors));
> PetscCall(DMDAVecGetArrayRead(cda, local, &coorslocal));
> PetscCall(DMDAGetCorners(cda, &mstart, &nstart, &pstart, &m, &n, &p));
> for (i = mstart; i < mstart + m; i++) {
> for (j = nstart; j < nstart + n; j++) {
> for (k = pstart; k < pstart + p; k++) {
> coors[k][j][i].x = coorslocal[k][j][i].x * F(0, 0);
> coors[k][j][i].y = coorslocal[k][j][i].y * F(1, 1);
> coors[k][j][i].z = coorslocal[k][j][i].z * F(2, 2);
> }
> }
> }
> PetscCall(DMDAVecRestoreArray(cda, global, &coors));
> PetscCall(DMDAVecRestoreArrayRead(cda, local, &coorslocal));
>
> PetscCall(DMGlobalToLocalBegin(cda, global, INSERT_VALUES, local));
> PetscCall(DMGlobalToLocalEnd(cda, global, INSERT_VALUES, local));
>
> PetscFunctionReturn(PETSC_SUCCESS);
> }
>
> On 17 Jan 2025, at 18:01, MIGUEL MOLINOS PEREZ <mmolinos at us.es> wrote:
>
> You are right!! Thank you again!
>
> Miguel
>
> On Jan 17, 2025, at 5:18 PM, Matthew Knepley <knepley at gmail.com> wrote:
>
> On Fri, Jan 17, 2025 at 10:49 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es>
> wrote:
>
>> Now the error is in the call to DMSwarmMigrate
>>
>
> You have almost certainly overwritten memory somewhere. Can you use
> vlagrind or Address Sanitizer?
>
> Thanks,
>
> Matt
>
>
>> Miguel
>>
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
>> probably memory access out of range
>> [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
>> [0]PETSC ERROR: or see https://urldefense.us/v3/__https://petsc.org/release/faq/*valgrind__;Iw!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpQ8EGYjN$ and
>> https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpcRZYRM5$
>> [0]PETSC ERROR: --------------------- Stack Frames
>> ------------------------------------
>> [0]PETSC ERROR: The line numbers in the error traceback are not always
>> exact.
>> [0]PETSC ERROR: #1 DMSwarmDataBucketGetSizes() at
>> /Users/migmolper/petsc/src/dm/impls/swarm/data_bucket.c:297
>> [0]PETSC ERROR: #2 DMSwarmMigrate_CellDMScatter() at
>> /Users/migmolper/petsc/src/dm/impls/swarm/swarm_migrate.c:201
>> [0]PETSC ERROR: #3 DMSwarmMigrate() at
>> /Users/migmolper/petsc/src/dm/impls/swarm/swarm.c:1349
>> [0]PETSC ERROR: #4 main() at
>> /Users/migmolper/DMD/driver-tasting-SOLERA.cpp:41
>> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
>>
>> On Jan 17, 2025, at 4:22 PM, Matthew Knepley <knepley at gmail.com> wrote:
>>
>> On Fri, Jan 17, 2025 at 10:08 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es>
>> wrote:
>>
>>> Thank you Matt, this the piece of code I use to change the coordinates
>>> of the DM obtained using:
>>>
>>
>> You do not need the call to DMSetCoordinates(). What happens when you
>> remove it?
>>
>> Thanks,
>>
>> Matt
>>
>>
>>>
>>> DMSwarmGetCellDM(Simulation.atomistic_data, &bounding_cell);
>>> DMGetApplicationContext(bounding_cell, &background_mesh);
>>>
>>> Thanks,
>>> Miguel
>>>
>>>
>>> /************************************************************************/
>>>
>>> PetscErrorCode Volumetric_Expansion(DM dm, const Eigen::Matrix3d& F) {
>>> PetscErrorCode ierr;
>>> Vec coordinates;
>>> PetscScalar* coordArray;
>>> PetscInt xs, ys, zs, xm, ym, zm, i, j, k;
>>> PetscInt dim, M, N, P;
>>>
>>> PetscFunctionBegin;
>>> // Get DMDA information
>>> ierr = DMDAGetInfo(dm, &dim, &M, &N, &P, NULL, NULL, NULL, NULL, NULL,
>>> NULL,
>>> NULL, NULL, NULL);
>>> CHKERRQ(ierr);
>>> ierr = DMDAGetCorners(dm, &xs, &ys, &zs, &xm, &ym, &zm);
>>> CHKERRQ(ierr);
>>>
>>> // Get the coordinates vector
>>> ierr = DMGetCoordinates(dm, &coordinates);
>>> CHKERRQ(ierr);
>>> ierr = VecGetArray(coordinates, &coordArray);
>>> CHKERRQ(ierr);
>>>
>>> // Update the coordinates based on the desired transformation
>>> for (k = zs; k < zs + zm; k++) {
>>> for (j = ys; j < ys + ym; j++) {
>>> for (i = xs; i < xs + xm; i++) {
>>> PetscInt idx =
>>> ((k * N + j) * M + i) * dim; // Index for the i, j, k point
>>> coordArray[idx] = coordArray[idx] * F(0,0); // Update x-coordinate
>>> coordArray[idx + 1] = coordArray[idx + 1] * F(1,1); // Update
>>> y-coordinate
>>> coordArray[idx + 2] = coordArray[idx + 2] * F(2,2); // Update
>>> z-coordinate
>>> }
>>> }
>>> }
>>>
>>> // Restore the coordinates vector
>>> ierr = VecRestoreArray(coordinates, &coordArray);
>>> CHKERRQ(ierr);
>>>
>>> // Set the updated coordinates back to the DMDA
>>> ierr = DMSetCoordinates(dm, coordinates);
>>> CHKERRQ(ierr);
>>>
>>> PetscFunctionReturn(0);
>>> }
>>>
>>>
>>> /************************************************************************/
>>>
>>> On 17 Jan 2025, at 16:00, Matthew Knepley <knepley at gmail.com> wrote:
>>>
>>> On Fri, Jan 17, 2025 at 9:45 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es>
>>> wrote:
>>>
>>>> I tried what you suggested, but still I got this error message. Maybe I
>>>> should use main release?
>>>>
>>>
>>> No. I suspect something is wrong with the way you are setting
>>> coordinates. Can you share the code?
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>>
>>>> Miguel
>>>>
>>>> [4]PETSC ERROR:
>>>> ------------------------------------------------------------------------
>>>> [4]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
>>>> probably memory access out of range
>>>> [4]PETSC ERROR: Try option -start_in_debugger or
>>>> -on_error_attach_debugger
>>>> [4]PETSC ERROR: or see https://urldefense.us/v3/__https://petsc.org/release/faq/*valgrind__;Iw!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpQ8EGYjN$ and
>>>> https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpcRZYRM5$
>>>> [4]PETSC ERROR: --------------------- Stack Frames
>>>> ------------------------------------
>>>> [4]PETSC ERROR: The line numbers in the error traceback are not always
>>>> exact.
>>>> [4]PETSC ERROR: #1 Pack_PetscReal_1_0() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/impls/basic/sfpack.c:373
>>>> [4]PETSC ERROR: #2 PetscSFLinkPackRootData_Private() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/impls/basic/sfpack.c:932
>>>> [4]PETSC ERROR: #3 PetscSFLinkPackRootData() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/impls/basic/sfpack.c:966
>>>> [4]PETSC ERROR: #4 PetscSFBcastBegin_Basic() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/impls/basic/sfbasic.c:357
>>>> [4]PETSC ERROR: #5 PetscSFBcastWithMemTypeBegin() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/interface/sf.c:1513
>>>> [4]PETSC ERROR: #6 VecScatterBegin_Internal() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/interface/vscat.c:70
>>>> [4]PETSC ERROR: #7 VecScatterBegin() at
>>>> /Users/migmolper/petsc/src/vec/is/sf/interface/vscat.c:1316
>>>> [4]PETSC ERROR: #8 DMGlobalToLocalBegin_DA() at
>>>> /Users/migmolper/petsc/src/dm/impls/da/dagtol.c:15
>>>> [4]PETSC ERROR: #9 DMGlobalToLocalBegin() at
>>>> /Users/migmolper/petsc/src/dm/interface/dm.c:2844
>>>> [4]PETSC ERROR: #10 DMGetCoordinatesLocalSetUp() at
>>>> /Users/migmolper/petsc/src/dm/interface/dmcoordinates.c:565
>>>> [4]PETSC ERROR: #11 DMGetCoordinatesLocal() at
>>>> /Users/migmolper/petsc/src/dm/interface/dmcoordinates.c:599
>>>> [4]PETSC ERROR: #12 _DMLocatePoints_DMDARegular_IS() at
>>>> /Users/migmolper/DMD/SOLERA/Atoms/Atom.cpp:531
>>>> [4]PETSC ERROR: #13 DMLocatePoints_DMDARegular() at
>>>> /Users/migmolper/DMD/SOLERA/Atoms/Atom.cpp:586
>>>> [4]PETSC ERROR: #14 DMLocatePoints() at
>>>> /Users/migmolper/petsc/src/dm/interface/dmcoordinates.c:1194
>>>> [4]PETSC ERROR: #15 DMSwarmMigrate_CellDMScatter() at
>>>> /Users/migmolper/petsc/src/dm/impls/swarm/swarm_migrate.c:219
>>>> [4]PETSC ERROR: #16 DMSwarmMigrate() at
>>>> /Users/migmolper/petsc/src/dm/impls/swarm/swarm.c:1349
>>>> [4]PETSC ERROR: #17 main() at
>>>> /Users/migmolper/DMD/driver-tasting-SOLERA.cpp:41
>>>>
>>>>
>>>>
>>>> On Jan 15, 2025, at 4:56 PM, MIGUEL MOLINOS PEREZ <mmolinos at us.es>
>>>> wrote:
>>>>
>>>> Thank you Matt for the useful info. I’ll try your idea.
>>>>
>>>> Miguel
>>>>
>>>> On 15 Jan 2025, at 16:48, Matthew Knepley <knepley at gmail.com> wrote:
>>>>
>>>> On Wed, Jan 15, 2025 at 10:41 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es>
>>>> wrote:
>>>>
>>>>> Thank you Matt.
>>>>>
>>>>> Yes, I am getting the "CellDM" from the DMSwarm.
>>>>>
>>>>> 1. I have recently overhauled this functionality because it was not
>>>>> flexible enough for the plasma simulation we do. Thus main and release work
>>>>> differently.
>>>>>
>>>>>
>>>>> Nice to hear that. Should I move to main?
>>>>>
>>>>
>>>> The changes allow you to have several cell DMs. I want to bin particles
>>>> in space, but also in velocity, and then in the tensor product of space and
>>>> velocity. Moreover, sometimes I want to use different Swarm fields as the
>>>> DM field for the solver. You can do all that with main now. If you just
>>>> need a single DM with the same DM fields, release is fine.
>>>>
>>>>
>>>>> 2. I assume you are using release
>>>>>
>>>>>
>>>>> You are correct.
>>>>>
>>>>> 3. In both main and release, if you change the coordinates of your
>>>>> CellDM mesh, you need to rebin the particles. The easiest way to do this is
>>>>> to call DMSwarmMigrate(sw, PETSC_FALSE).
>>>>>
>>>>>
>>>>> What do you mean by rebin?
>>>>>
>>>>
>>>> When you provide the cell DM, Swrm makes a "sort context" that bins the
>>>> particles into DM cells. If you change the coordinates, this binning will
>>>> change, so you need it to "rebin" or recreate the sort context.
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>>
>>>>> Miguel
>>>>>
>>>>>
>>>>> Thanks,
>>>>>
>>>>> Matt
>>>>>
>>>>>
>>>>>> Best,
>>>>>> Miguel
>>>>>>
>>>>> --
>>>>> 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
>>>>>
>>>>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$
>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> 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
>>>>
>>>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$
>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
>>>>
>>>>
>>>>
>>>>
>>>
>>> --
>>> 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
>>>
>>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$
>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
>>>
>>>
>>>
>>
>> --
>> 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
>>
>> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
>>
>>
>>
>
> --
> 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
>
> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
>
>
>
>
--
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
https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpSDqciPi$ <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d_kL76N30FqcKUQr7kFkiAg2Rt5XrsdTHF_opDlqKEu3B89-dlXJyQgH64brfGhgW5PJg0jmxzkTpUuriFiH$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250124/3f8e2ac7/attachment-0001.html>
More information about the petsc-users
mailing list