<div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote">On Thu, Nov 14, 2013 at 1:05 PM, Bishesh Khanal <span dir="ltr"><<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><br><div class="gmail_extra"><br><br><div class="gmail_quote"><div class="im">On Wed, Nov 13, 2013 at 10:06 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:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div>Bishesh Khanal <<a href="mailto:bisheshkh@gmail.com" target="_blank">bisheshkh@gmail.com</a>> writes:<br>
<br>
> Within A, for now, I can consider mu to be constant, although later if<br>
> possible it can be a variable even a tensor to describe anisotropy. But to<br>
> start with I want put this a constant.<br>
> The original equations start with mu (grad(u) + grad(u)^T) but then<br>
> simplifications occur due to div(u) = f2<br>
<br>
</div>Rework that step in case of variable mu.<br></blockquote><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div><br>
> I'm mostly interested in the phenomenon in A with my model, here B is the<br>
> extension of the very irregular domain of A to get a cuboid. Here, in B I<br>
> release the div(u) = f2 constraint and just put a regularisation to<br>
> penalize large deformation. What is of importance here is to compensate the<br>
> net volume expansion in domain A by corresponding contraction in domain B<br>
> so that the boundaries of the cuboid do not move. It does not particularly<br>
> represent any physics except probably that it gives me a velocity field<br>
> having a certain divergence field that penalizes big deformations.<br>
<br>
</div>Okay, sounds like it's already an artificial equation, so you should be<br>
able to leave in a normal equation for p, with a big mass matrix on the<br>
diagonal,<br>
<div><br>
div(mu(grad(u))) - grad(p) = f1<br>
</div>div(u) - c(x) p = f2<br>
<br>
c(x) = 0 in domain A and c(x) is large (the inverse of the second Lamé<br>
parameter) in domain B.<br>
<div><br></div></blockquote></div></div></div></div></blockquote><div><br></div><div>I tried to modify my existing code to incorporate c(x) instead of 0 as you suggested. When c(x) is zero everywhere, the solver converges and gives me the result I expect, as before. However, when c(x) is non-zero in domain B, the solver does not seem to converge. The ouput for both cases are given below. When c(x) was zero everywhere, I had set the constant pressure null space to the system matrix and the schur complement matrix since I am solving the system with pctype fieldsplit and schur. However, when c(x) is non-zero in domain B (I put c(x) = 1 in B for test), we do not have the constant null space for pressure anymore; so for this case I do not set the constant null space. Can you comment on what could be possibly going wrong ? <br>
</div><div>1. Output for the test case: c(x) = 0 in A, c(x) = 1 in B; mu = 1 everywhere.<br>options given: <br>-pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_dm_splits 0 -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_0_pc_factor_mat_solver_package superlu_dist -fieldsplit_0_ksp_monitor -fieldsplit_0_ksp_converged_reason -fieldsplit_0_ksp_max_it 100 -fieldsplit_1_ksp_monitor -fieldsplit_1_ksp_converged_reason -ksp_converged_reason<br>
<br>output:<br><br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 2.321914527334e+01 <br> 1 KSP Residual norm 5.204444045635e-14 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
Residual norms for fieldsplit_1_ solve.<br> 0 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.993511593752e+00 <br> 1 KSP Residual norm 4.050802627628e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 1 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 6.284336607097e-01 <br> 1 KSP Residual norm 2.121905965324e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 2 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.182615416676e+00 <br> 1 KSP Residual norm 6.481437789588e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 3 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.470958139659e+00 <br> 1 KSP Residual norm 9.105094653121e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 4 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.802783367675e+00 <br> 1 KSP Residual norm 1.307376269273e-14 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 5 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 9.352429007750e-01 <br> 1 KSP Residual norm 5.214147133013e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 6 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 9.088208637857e-01 <br> 1 KSP Residual norm 5.202918034148e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 7 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 8.607808656970e-01 <br> 1 KSP Residual norm 5.212231778254e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 8 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.228224544246e+00 <br> 1 KSP Residual norm 7.032688183318e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 9 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.238989110105e+00 <br> 1 KSP Residual norm 7.308937417623e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 10 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 9.315232322259e-01 <br> 1 KSP Residual norm 4.434343770858e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 11 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.070108662141e+00 <br> 1 KSP Residual norm 5.822139925590e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 12 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.066025697964e+00 <br> 1 KSP Residual norm 6.055553792441e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 13 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 7.981752531989e-01 <br> 1 KSP Residual norm 3.750650135599e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 14 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 7.225821573201e-01 <br> 1 KSP Residual norm 3.043640976617e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 15 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 8.223502060699e-01 <br> 1 KSP Residual norm 4.212671738422e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 16 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 8.322451979979e-01 <br> 1 KSP Residual norm 2.999201599777e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 17 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.052268760248e+00 <br> 1 KSP Residual norm 5.911203456590e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 18 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.082754048106e+00 <br> 1 KSP Residual norm 5.992300533204e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br> 19 KSP Residual norm 2.890735677419e+14 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.061699839350e+00 <br>^C[mpiexec@edwards] Sending Ctrl-C to processes as requested<br>
</div><div>and it keeps on going like this.<br><br></div><div>The system is converging when I set incompressibility constraint everywhere. Here is a sample output for this case: c(x) = 0 everywhere, mu = 1 everywhere.<br>
</div><div>Options used:<br>-pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_dm_splits 0 -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 3 -fieldsplit_0_pc_type lu -fieldsplit_0_pc_factor_mat_solver_package superlu_dist -fieldsplit_0_ksp_monitor -fieldsplit_0_ksp_converged_reason -fieldsplit_0_ksp_max_it 100 -fieldsplit_1_ksp_monitor -fieldsplit_1_ksp_converged_reason -ksp_converged_reason<br>
<br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 0.000000000000e+00 <br> Linear solve converged due to CONVERGED_ATOL iterations 0<br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 2.321914527334e+01 <br>
1 KSP Residual norm 5.204444045635e-14 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> Residual norms for fieldsplit_1_ solve.<br> 0 KSP Residual norm 9.964395910568e+14 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.573745221888e+00 <br> 1 KSP Residual norm 3.612233230950e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 1 KSP Residual norm 2.667315756754e+12 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.053280683677e+00 <br> 1 KSP Residual norm 5.770507184021e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 2 KSP Residual norm 3.176176643790e+11 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.016103066894e+00 <br> 1 KSP Residual norm 4.558255925705e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 3 KSP Residual norm 1.046851012581e+11 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.095000885520e+00 <br> 1 KSP Residual norm 6.356842242710e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 4 KSP Residual norm 4.078940666596e+10 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.244953433655e+00 <br> 1 KSP Residual norm 6.745548989971e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 5 KSP Residual norm 1.774320878731e+10 <br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 1.176603450320e+00 <br> 1 KSP Residual norm 5.385862671655e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 6 KSP Residual norm 8.586271019948e+09 <br> Linear solve converged due to CONVERGED_RTOL iterations 6<br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.162170277427e+01 <br> 1 KSP Residual norm 2.912041048579e-14 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> Residual norms for fieldsplit_0_ solve.<br>
0 KSP Residual norm 9.290138250365e-01 <br> 1 KSP Residual norm 2.116243919788e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> Residual norms for fieldsplit_1_ solve.<br> 0 KSP Residual norm 3.986822705891e+13 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.573745244651e+00 <br> 1 KSP Residual norm 3.715168656229e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 1 KSP Residual norm 1.067205796344e+11 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.053498210605e+00 <br> 1 KSP Residual norm 4.526737475809e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 2 KSP Residual norm 1.270173568462e+10 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.017546903473e+00 <br> 1 KSP Residual norm 4.593425818353e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 3 KSP Residual norm 4.160360750454e+09 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.097951946288e+00 <br> 1 KSP Residual norm 6.461246234077e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 4 KSP Residual norm 1.609903325354e+09 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.250441943430e+00 <br> 1 KSP Residual norm 7.157658262436e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 5 KSP Residual norm 6.919994735907e+08 <br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.169858508137e+00 <br> 1 KSP Residual norm 5.046051842799e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> 6 KSP Residual norm 3.217042876248e+08 <br>
Linear solve converged due to CONVERGED_RTOL iterations 6<br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 4.649924240988e-01 <br> 1 KSP Residual norm 1.068254961936e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 7.194107998977e-04 <br> 1 KSP Residual norm 1.744183339526e-18 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br> Residual norms for fieldsplit_1_ solve.<br>
0 KSP Residual norm 4.448127329286e+12 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 7.865130501209e-01 <br> 1 KSP Residual norm 5.153218330579e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
1 KSP Residual norm 2.638913823938e+12 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.283369467124e+00 <br> 1 KSP Residual norm 7.504488207954e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
2 KSP Residual norm 1.426102610607e+12 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.129662815202e+00 <br> 1 KSP Residual norm 5.995004855070e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
3 KSP Residual norm 8.025360882010e+11 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.424860285739e+00 <br> 1 KSP Residual norm 8.619123051910e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
4 KSP Residual norm 5.286840519981e+11 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.289665727094e+00 <br> 1 KSP Residual norm 7.046395176709e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
5 KSP Residual norm 3.634745308165e+11 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.364708363745e+00 <br> 1 KSP Residual norm 7.402799565707e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
6 KSP Residual norm 1.720693722234e+11 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.102971050925e+00 <br> 1 KSP Residual norm 6.835965584112e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
7 KSP Residual norm 7.870496175081e+10 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.345118128097e+00 <br> 1 KSP Residual norm 7.734345056964e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
8 KSP Residual norm 2.921176674626e+10 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.291212792694e+00 <br> 1 KSP Residual norm 7.780332549058e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
9 KSP Residual norm 1.148381817318e+10 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.073046392930e+00 <br> 1 KSP Residual norm 6.464513846538e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
10 KSP Residual norm 5.253676669489e+09 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.016140779574e+00 <br> 1 KSP Residual norm 5.990033309135e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
11 KSP Residual norm 2.216200526999e+09 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.037961855022e+00 <br> 1 KSP Residual norm 6.208761778424e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
12 KSP Residual norm 9.013987504565e+08 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.033050143801e+00 <br> 1 KSP Residual norm 5.853852773652e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
13 KSP Residual norm 3.751638125875e+08 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 9.581056747534e-01 <br> 1 KSP Residual norm 4.028342105749e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
14 KSP Residual norm 1.372411728307e+08 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 1.117841059532e+00 <br> 1 KSP Residual norm 4.848302320534e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
15 KSP Residual norm 5.353181748969e+07 <br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 8.645767106479e-01 <br> 1 KSP Residual norm 4.186764891168e-15 <br> Linear solve converged due to CONVERGED_RTOL iterations 1<br>
16 KSP Residual norm 2.421357309016e+07 <br> Linear solve converged due to CONVERGED_RTOL iterations 16<br> Residual norms for fieldsplit_0_ solve.<br> 0 KSP Residual norm 6.302684114623e-01 <br> 1 KSP Residual norm 2.299313600630e-15 <br>
Linear solve converged due to CONVERGED_RTOL iterations 1<br>Linear solve converged due to CONVERGED_RTOL iterations 2<br></div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote"><div class="im"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div></div></blockquote>
</div><div>Thanks, this looks quite reasonable. I'll try to experiment with it.<br></div><div class="im"><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<div>
> I do not know much about FEM. But some of the reasons why I have avoided it<br>
> in this particular problem are: (Please correct me on any of the following<br>
> points if they are wrong)<br>
> 1. The inputs f1 and f2 are 3D images (in average of size 200^3) that come<br>
> from other image processing pipeline; it's important that I constrain u at<br>
> each voxel for div(u) = f2 in domain A. I am trying to avoid having to get<br>
> the meshing from the 3D image(with very detailed structures), then go back<br>
> to the image from the obtained u again because I have to use the obtained u<br>
> to warp the image, transport other parameters again with u in the image<br>
> space and again obtain new f1 and f2 images. Then iterate this few times.<br>
<br>
</div>Okay, there's nothing wrong with that.<br>
</blockquote></div></div>Thanks for the confirmation.<br></div><div class="gmail_extra"><br></div></div>
</blockquote></div><br></div></div>