[petsc-users] solving stokes-like equation in 3D staggered grid in irregular domain with petscsection

[Bishesh Khanal]

On Tue, Nov 12, 2013 at 10:55 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:

> Bishesh Khanal <bisheshkh at gmail.com> writes:
>
> > Dear all,
> > I have an implementation of stokes flow equation in cuboid using a
> > staggered grid with DMDA. Now I want to solve a problem in a slightly
> > different (but more difficult) situation where I think PetscSection might
> > be useful although I have not used PetscSection that much and I
> understand
> > that some of its aspects is still under development in petsc.
> > I'm ok to work with the dev version of petsc if what I want to do using
> > petscsection is plausible.
> >
> > Here is my problem description:
> >
> > I have a 3D binary image that partitions a cuboid into domain A and
> domain
> > B by tagging each cell as 0 or 1. The shape of A and therefore B can be
> > fairly arbitrary. I want to solve the following sets of equations:
> >
> > A. In A:
> > div(mu(grad(u))) - grad(p) = f1   , where f1 depends on the position x.
>
> Is mu constant or variable?  Should this be mu (grad(u) + grad(u)^T)/2 ?
>

Within A, for now, I can consider mu to be constant, although later if
possible it can be a variable even a tensor to describe anisotropy. But to
start with I want put this a constant.
The original equations start with mu (grad(u) + grad(u)^T) but then
simplifications occur due to div(u) = f2 where f2 is known and is an input
to the system. Some terms involving f2 go to the right hand side to form f1.


>
> > div(u)                                   = f2   , where f2 is non-zero
> and
> > depends on the position x.
>
> Does it really only depend on position or does it depend on state
> variables (u,p)?
>

f2 comes from my other models, it is an output of an image processing
pipeline, thus is a 3D image. So it does vary depending on position, and is
independent of (u,p) for this time independent sets of equations.


>
> > ---------------------------------------
> > B. In B:
> > div(mu(grad(u)) = f3,    where f3 depends on the position x.
>
> What physics does this represent?
>

I'm mostly interested in the phenomenon in A with my model, here B is the
extension of the very irregular domain of A to get a cuboid. Here, in B I
release the div(u) = f2 constraint and just put a regularisation to
penalize large deformation. What is of importance here is to compensate the
net volume expansion in domain A by corresponding contraction in domain B
so that the boundaries of the cuboid do not move. It does not particularly
represent any physics except probably that it gives me a velocity field
having a certain divergence field that penalizes big deformations.


> > -----------------------------------------
> > u = constant on the cuboid boundaries. (Dirichlet boundaries)
> > p = constant on domain B (if needed as boundary condition for domain A)
> > -------------------------------------------
> >
> > Now, I was thinking of creating a single matrix M that contains the
> > operators for both A and B domains, and solve a linear system MX = R,
> where
> > X and R are the solution and Rhs vectors of roughly the sizes 4*nA + 3*nB
> > where,   nA => no. of cells in A, nB => no. of cells in B.
> >
> > Some questions:
> >
> > 1. Would PetscSection be useful in this case due to the arbitrary shape
> of
> > domain A and because of different number of dofs in domain A and domain
> B ?
> > Is there a simple example that uses petscsection for the staggered grid
> > case ? I would like to be able to implement staggered grid without the
> use
> > of ghost/fictitious cells.
>
> Just put degrees of freedom on faces and cell centers when creating the
> section.
>
>
I've not understood petscsection well. I tried reading up in the manual but
it does not look easy to follow with all the terms I'm unfamiliar with. Is
it possible to have an example for the staggered grid case ? Or a short
illustration with code snippets/idea on how should I start with it ?


> > 2. Is it a good idea to create a single matrix corresponding to both the
> > domains A and B (particularly in view of the difficulty it might bring in
> > solving with some preconditioners using schur fieldsplit) ?
>
> I would start with a single matrix.  Use MatSetValuesLocal() in assembly
> so that you can easily switch to a different format.
>
> > 3. I'm trying to put both domains in a single matrix to avoid the
> > difficulty I would have if I want to consider only the domain A. In this
> > case I would need a traction free boundary condition on the irregular
> > boundary of domain A, and it seems a bit too challenging for me to
> > incorporate it with the staggered grid. If there is an idea to implement
> > this and if you think this could be more suitable than the approach in 2
> > above, I would like to learn about that too!
>
> Complexity of implementing boundary conditions on staggered grids is one
> reason some people turn to other discretization technology, such as
> finite elements.
>

I do not know much about FEM. But some of the reasons why I have avoided it
in this particular problem are:  (Please correct me on any of the following
points if they are wrong)
1. The inputs f1 and f2 are 3D images (in average of size 200^3) that come
from other image processing pipeline; it's important that I constrain u at
each voxel for div(u) = f2 in domain A. I am trying to avoid having to get
the meshing from the 3D image(with very detailed structures), then go back
to the image from the obtained u again because I have to use the obtained u
to warp the image, transport other parameters again with u in the image
space and again obtain new f1 and f2 images. Then iterate this few times.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131113/58915766/attachment-0001.html>