Sorry for the delay. The fix is in the git branch barry/2024-01-09/fix-mirror-dmda-3d/release see also https://gitlab.com/petsc/petsc/-/merge_requests/7175 <ttps://gitlab.com/petsc/petsc/-/merge_requests/7175>
Barry > On Jan 8, 2024, at 3:44 PM, Gourav Kumbhojkar <[email protected]> > wrote: > > You are right. Attaching the code that I’m using to test this. > The output matrix is saved in separate ascii files. > You can use “make noflux” to compile the code. > > Gourav > > From: Barry Smith <[email protected] <mailto:[email protected]>> > Date: Saturday, January 6, 2024 at 7:08 PM > To: Gourav Kumbhojkar <[email protected] > <mailto:[email protected]>> > Cc: [email protected] <mailto:[email protected]> > <[email protected] <mailto:[email protected]>> > Subject: Re: [petsc-users] Neumann Boundary Condition with DMDACreate3D > > > If the mirror code for star stencil is just wrong in 3d we should simply > fix it. Not use some other approach. Can you attach code that tries to do > what you need for both 2d (that results in a matrix you are happy with) and > 3d (that results in a matrix that you are not happy with). > > Barry > > > > > > On Jan 6, 2024, at 7:30 PM, Gourav Kumbhojkar <[email protected] > <mailto:[email protected]>> wrote: > > Thank you, Barry. Sorry for the late response. > > Yes, I was referring to the same text. I am using a star stencil. However, I > don’t think the mirror condition is implemented for star stencil either. > > TLDR version of the whole message typed below – > I think DM_BOUNDARY_GHOSTED is not implemented correctly in 3D. It appears > that ghost nodes are mirrored with boundary nodes themselves. They should > mirror with the nodes next to boundary. > > Long version - > Here’s what I’m trying to do – > > Step 1 - Create a 3D DM > ierr = DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_MIRROR, DM_BOUNDARY_MIRROR, > DM_BOUNDARY_MIRROR, DMDA_STENCIL_STAR, num_pts, num_pts, num_pts, > PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, NULL, &da); > CHKERRQ(ierr); > Note - num_pts = 4 in my code. > > Step 2 – Create a Matrix from DM ( a FDM stiffness matrix) > DMCreateMatrix(da, &K); > globalKMat(K, info); > > “globalKMat” is a user-defined function. Here’s a snippet from this function: > for (int i = info.xs; i < (info.xs + info.xm); i++){ > for(int j = info.ys; j < (info.ys + info.ym); j++){ > for (int k = info.zs; k < (info.zs + info.zm); k++){ > ncols = 0; > row.i = i; row.j = j; row.k = k; > > col[0].i = i; col[0].j = j; col[0].k = k; > vals[ncols++] = -6.; //ncols=1 > > col[ncols].i = i-1; col[ncols].j = j; col[ncols].k = k; > vals[ncols++] = 1.;//ncols=2 > > There are total 7 “ncols”. Other than the first one all ncols have value 1 > (first one is set to -6). As you can see, this step is to only build the FDM > stiffness matrix. I use “ADD_VALUES” at the end in the above function. > > Step 3 – View the stiffness matrix to check the values. I use MatView for > this. > > Here are the results – > 3D DM (showing first few rows of the stiffness matrix here, the original > matrix is 64x64)- > Mat Object: 1 MPI processes > type: seqaij > row 0: (0, -3.) (1, 1.) (4, 1.) (16, 1.) > row 1: (0, 1.) (1, -4.) (2, 1.) (5, 1.) (17, 1.) > row 2: (1, 1.) (2, -4.) (3, 1.) (6, 1.) (18, 1.) > row 3: (2, 1.) (3, -3.) (7, 1.) (19, 1.) > row 4: (0, 1.) (4, -4.) (5, 1.) (8, 1.) (20, 1.) > row 5: (1, 1.) (4, 1.) (5, -5.) (6, 1.) (9, 1.) (21, 1.) > > Repeat the same steps for a 2D DM to show the difference (the entire matrix > is now 16x16) > Mat Object: 1 MPI processes > type: seqaij > row 0: (0, -4.) (1, 2.) (4, 2.) > row 1: (0, 1.) (1, -4.) (2, 1.) (5, 2.) > row 2: (1, 1.) (2, -4.) (3, 1.) (6, 2.) > row 3: (2, 2.) (3, -4.) (7, 2.) > row 4: (0, 1.) (4, -4.) (5, 2.) (8, 1.) > row 5: (1, 1.) (4, 1.) (5, -4.) (6, 1.) (9, 1.) > > I suspect that when using “DM_BOUNDARY_MIRROR” in 3D, the ghost node value is > added to the boundary node itself, which would explain why row 0 of the > stiffness matrix has -3 instead of -6. In principle the ghost node value > should be mirrored with the node next to boundary. > Clearly, there’s no issue with the 2D implementation of the mirror boundary. > The row 0 values are -4, 2, and 2 as expected. > > Let me know if I should give any other information about this. I also thought > about using DM_BOUNDARY_GHOSTED and implement the mirror boundary in 3D from > scratch but I would really appreciate some resources on how to do that. > > Thank you. > > Gourav > > > From: Barry Smith <[email protected] <mailto:[email protected]>> > Date: Thursday, January 4, 2024 at 12:24 PM > To: Gourav Kumbhojkar <[email protected] > <mailto:[email protected]>> > Cc: [email protected] <mailto:[email protected]> > <[email protected] <mailto:[email protected]>> > Subject: Re: [petsc-users] Neumann Boundary Condition with DMDACreate3D > > > Are you referring to the text? > > . `DM_BOUNDARY_MIRROR` - the ghost value is the same as the value 1 grid > point in; that is, the 0th grid point in the real mesh acts like a mirror to > define > the ghost point value; not yet implemented for 3d > > > Looking at the code for DMSetUp_DA_3D() I see > > PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && > by != DM_BOUNDARY_MIRROR && bz != DM_BOUNDARY_MIRROR), > PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box > stencil"); > > which seems (to me) to indicate the mirroring is not done for box stencils > but should work for star stencils. > > Are you using a star stencil or a box stencil? > > I believe the code is not complete for box stencil because the code to > determine the location of the "mirrored point" for extra "box points" is > messy in 3d and no one wrote it. You can compare DMSetUp_DA_2D() and > DMSetUp_DA_3D() to see what is missing and see if you can determine how to > add it for 3d. > > Barry > > > > On Jan 4, 2024, at 1:08 PM, Gourav Kumbhojkar <[email protected] > <mailto:[email protected]>> wrote: > > Hi, > > I am trying to implement a No-flux boundary condition for a 3D domain. I > previously modeled a no flux boundary in 2D domain using DMDACreate2D and > “PETSC_BOUNDARY_MIRROR” which worked well. > However, the manual pages say that the Mirror Boundary is not supported for > 3D. > Could you please point me to the right resources to implement no flux > boundary condition in 3D domains? > > Regards, > Gourav K. > > <noflux_check.c><makefile>
