Solving Poisson eqn wth Dirichlet boundaries inside the domain

[Mani Chandra]

Hello,

I'm trying to solve the Poisson eqn with Dirichlet boundaries inside the domain. I have a structured DM grid coupled to SNES as follows:

da = PETSc.DMDA().create([N1, N2],
                         stencil_width = N_ghost,
                         boundary_type = ('periodic',
                                          'periodic'
                                         ),
                         stencil_type = 1,
                         dof          = 1,
                         comm         = comm
                       )

glob_residual  = da.createGlobalVec()
glob_phi = da.createLocalVec()

snes = PETSc.SNES().create()
pde  = poisson_eqn()
snes.setFunction(pde.compute_residual, glob_residual)

snes.setDM(da) # Allows for use of graph coloring during Jacobian assembly
snes.setFromOptions()

snes.solve(None, glob_phi)

Now to specify the Dirichlet boundaries on a surface inside the domain, I'm trying the following inside the residual assembly function (pde.compute_residual):

phi[interior_indices] = 1. # Dirichlet BCs in the interior of the domain
residual[interior_indices] = 0 # Remove contribution to the residual norm from this surface

However, this approach doesn't seem to be working with SNES not converging:

    snes.solve(None, glob_phi)
  File "PETSc/SNES.pyx", line 520, in petsc4py.PETSc.SNES.solve (src/petsc4py.PETSc.c:167948)
petsc4py.PETSc.Error: error code 91
[0] SNESSolve() line 4017 in /home/mani/Downloads/petsc-3.7.5/src/snes/interface/snes.c
[0]   
[0] SNESSolve has not converged

Is there a way to correctly solve the problem of imposing Dirichlet boundaries in the interior of the domain?

Thanks,

Mani

[Lisandro Dalcin]

Look at demo/bratu3d/bratu3d.py, go to comment line # boundary points and the line f[i, j, k] = x[i, j, k] - 0, this example uses homogeneous Dirichlet BCs, otherwise you have to use a BC value other than 0.

PS: I would suggest to ask these kind of questions in the petsc-users mailing list.

[Mani Chandra]

I think there's some confusion here. I was asking about how to impose boundary conditions inside the domain. For example, to apply a fixed potential onto the plates of a capacitor, with the capacitor well inside the computational domain. I need to be able to remove the points in the interior where I'm applying the boundary condition from the solver. Is that possible using the petsc4py interface to petsc?

The example demo/bratu3d/bratu3d.py shows how to impose boundary conditions in the boundary points, not in the interior of the domain.

[Lisandro Dalcin]

In this particular case, whether it is possible to do what you want is not about petsc4py and its interface. It is about how you implement these BCs in the context of finite differences. DMDA is a data structure to hande a distributed structure grid and related core algorithms, no less, no more. The filling of residual vectors is up to the user. YOU have to decide how and where and what to put in every entry of the residual vector (and Jacobian matrix).

All that being said, in finite differences, imposing Dirichlet boundary conditions inside the domain is not different at all than doing it on the boundary. You just need to figure which are the global indices of these internal points, and put there the boundary condition equation.