[petsc-users] Update DMDA attached to DMSWARM

Matthew Knepley knepley at gmail.com
Fri Jan 24 09:50:24 CST 2025


On Fri, Jan 24, 2025 at 10:36 AM MIGUEL MOLINOS PEREZ <mmolinos at us.es>
wrote:

> Thanks Matt, I tried that too, and the problem remains. The box is updated
> only if I set no periodic bcc.
>

What do you mean by "The box is updated"? I am trying to understand how you
test things. Clearly the coordinates are updated,
even in the periodic case. Thus, I need to understand the test. Once we do
that, we can work backwards to the first broken thing.

  Thanks,

      Matt


> Miguel
>
> On 24 Jan 2025, at 14:20, Matthew Knepley <knepley at gmail.com> wrote:
>
> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupuIaVdEg$  and
>>> https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupsTWdVgG$ 
>>> [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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupuIaVdEg$  and
>>>>> https://urldefense.us/v3/__https://petsc.org/release/faq/__;!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupsTWdVgG$ 
>>>>> [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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
>>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>>>>>>
>>>>>>
>>>>>>
>>>>>
>>>>> --
>>>>> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
>>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>>>>>
>>>>>
>>>>>
>>>>>
>>>>
>>>> --
>>>> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
>>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>>>>
>>>>
>>>>
>>>
>>> --
>>> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
>>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>>>
>>>
>>>
>>
>> --
>> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
>> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>>
>>
>>
>>
>
> --
> 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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$ 
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
>
>
>

-- 
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!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSuphgzDu4D$  <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bl1IKWYD5vUTh2-bmOUYT72lASrjcNt_e-FTAgWDrFKIB9bal_DAR9Nx9MMl8oqM4OEMlkElerSupiN9pla9$ >
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250124/d07db27a/attachment-0001.html>


More information about the petsc-users mailing list