Ok, I now understand what you are reporting and am working on a fix.
From the manual pages I see
DMBoundaryType - Describes the choice for the filling of ghost cells on
physical domain boundaries.
Values:
+ `DM_BOUNDARY_NONE` - no ghost nodes
. `DM_BOUNDARY_GHOSTED` - ghost vertices/cells exist but aren't filled; you can
put values into them and then apply a stencil that uses those ghost locations
. `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
....
Developer Note:
Should `DM_BOUNDARY_MIRROR` have the same meaning with DMDA_Q0, that is a
staggered grid? In that case should the ghost point have the same value
as the 0th grid point where the physical boundary serves as the mirror?
I assume you are working with a "vertex valued" stencil as opposed to the
cell-centered values where mirror possibly has a different interpretation?
> 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>
