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 –
1. 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.)
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]>
Date: Thursday, January 4, 2024 at 12:24 PM
To: Gourav Kumbhojkar <[email protected]>
Cc: [email protected] <[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]>
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.
