<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Tue, Nov 12, 2013 at 10:55 PM, Jed Brown <span dir="ltr"><<a href="mailto:jedbrown@mcs.anl.gov" target="_blank">jedbrown@mcs.anl.gov</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im">Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com">bisheshkh@gmail.com</a>> writes:<br>
<br>
> Dear all,<br>
> I have an implementation of stokes flow equation in cuboid using a<br>
> staggered grid with DMDA. Now I want to solve a problem in a slightly<br>
> different (but more difficult) situation where I think PetscSection might<br>
> be useful although I have not used PetscSection that much and I understand<br>
> that some of its aspects is still under development in petsc.<br>
> I'm ok to work with the dev version of petsc if what I want to do using<br>
> petscsection is plausible.<br>
><br>
> Here is my problem description:<br>
><br>
> I have a 3D binary image that partitions a cuboid into domain A and domain<br>
> B by tagging each cell as 0 or 1. The shape of A and therefore B can be<br>
> fairly arbitrary. I want to solve the following sets of equations:<br>
><br>
> A. In A:<br>
> div(mu(grad(u))) - grad(p) = f1 , where f1 depends on the position x.<br>
<br>
</div>Is mu constant or variable? Should this be mu (grad(u) + grad(u)^T)/2 ?<br></blockquote><div><br></div><div>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.<br>
</div><div>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.<br>
</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="im"><br>
> div(u) = f2 , where f2 is non-zero and<br>
> depends on the position x.<br>
<br>
</div>Does it really only depend on position or does it depend on state<br>
variables (u,p)?<br></blockquote><div><br></div><div>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. <br>
</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="im"><br>
> ---------------------------------------<br>
> B. In B:<br>
> div(mu(grad(u)) = f3, where f3 depends on the position x.<br>
<br>
</div>What physics does this represent?<br></blockquote><div><br></div><div>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. <br>
<br></div><div></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div class="im"><br>
> -----------------------------------------<br>
> u = constant on the cuboid boundaries. (Dirichlet boundaries)<br>
> p = constant on domain B (if needed as boundary condition for domain A)<br>
> -------------------------------------------<br>
><br>
> Now, I was thinking of creating a single matrix M that contains the<br>
> operators for both A and B domains, and solve a linear system MX = R, where<br>
> X and R are the solution and Rhs vectors of roughly the sizes 4*nA + 3*nB<br>
> where, nA => no. of cells in A, nB => no. of cells in B.<br>
><br>
> Some questions:<br>
><br>
> 1. Would PetscSection be useful in this case due to the arbitrary shape of<br>
> domain A and because of different number of dofs in domain A and domain B ?<br>
> Is there a simple example that uses petscsection for the staggered grid<br>
> case ? I would like to be able to implement staggered grid without the use<br>
> of ghost/fictitious cells.<br>
<br>
</div>Just put degrees of freedom on faces and cell centers when creating the<br>
section.<br>
<div class="im"><br></div></blockquote><div><br></div><div>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 ?<br>
</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div class="im">
> 2. Is it a good idea to create a single matrix corresponding to both the<br>
> domains A and B (particularly in view of the difficulty it might bring in<br>
> solving with some preconditioners using schur fieldsplit) ?<br>
<br>
</div>I would start with a single matrix. Use MatSetValuesLocal() in assembly<br>
so that you can easily switch to a different format.<br>
<div class="im"><br>
> 3. I'm trying to put both domains in a single matrix to avoid the<br>
> difficulty I would have if I want to consider only the domain A. In this<br>
> case I would need a traction free boundary condition on the irregular<br>
> boundary of domain A, and it seems a bit too challenging for me to<br>
> incorporate it with the staggered grid. If there is an idea to implement<br>
> this and if you think this could be more suitable than the approach in 2<br>
> above, I would like to learn about that too!<br>
<br>
</div>Complexity of implementing boundary conditions on staggered grids is one<br>
reason some people turn to other discretization technology, such as<br>
finite elements.<br></blockquote><div> </div></div>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)<br>
</div><div class="gmail_extra">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.<br>
</div><br></div>