On Thu, Nov 14, 2013 at 1:05 PM, Bishesh Khanal <[email protected]> wrote:
> > > > On Wed, Nov 13, 2013 at 10:06 PM, Jed Brown <[email protected]> wrote: > >> Bishesh Khanal <[email protected]> writes: >> >> > 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 >> >> Rework that step in case of variable mu. >> > >> > 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. >> >> Okay, sounds like it's already an artificial equation, so you should be >> able to leave in a normal equation for p, with a big mass matrix on the >> diagonal, >> >> div(mu(grad(u))) - grad(p) = f1 >> div(u) - c(x) p = f2 >> >> c(x) = 0 in domain A and c(x) is large (the inverse of the second Lamé >> parameter) in domain B. >> >> 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 ? 1. Output for the test case: c(x) = 0 in A, c(x) = 1 in B; mu = 1 everywhere. options given: -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 output: Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 2.321914527334e+01 1 KSP Residual norm 5.204444045635e-14 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.993511593752e+00 1 KSP Residual norm 4.050802627628e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 1 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 6.284336607097e-01 1 KSP Residual norm 2.121905965324e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 2 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.182615416676e+00 1 KSP Residual norm 6.481437789588e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 3 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.470958139659e+00 1 KSP Residual norm 9.105094653121e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 4 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.802783367675e+00 1 KSP Residual norm 1.307376269273e-14 Linear solve converged due to CONVERGED_RTOL iterations 1 5 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 9.352429007750e-01 1 KSP Residual norm 5.214147133013e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 6 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 9.088208637857e-01 1 KSP Residual norm 5.202918034148e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 7 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 8.607808656970e-01 1 KSP Residual norm 5.212231778254e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 8 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.228224544246e+00 1 KSP Residual norm 7.032688183318e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 9 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.238989110105e+00 1 KSP Residual norm 7.308937417623e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 10 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 9.315232322259e-01 1 KSP Residual norm 4.434343770858e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 11 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.070108662141e+00 1 KSP Residual norm 5.822139925590e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 12 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.066025697964e+00 1 KSP Residual norm 6.055553792441e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 13 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 7.981752531989e-01 1 KSP Residual norm 3.750650135599e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 14 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 7.225821573201e-01 1 KSP Residual norm 3.043640976617e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 15 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 8.223502060699e-01 1 KSP Residual norm 4.212671738422e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 16 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 8.322451979979e-01 1 KSP Residual norm 2.999201599777e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 17 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.052268760248e+00 1 KSP Residual norm 5.911203456590e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 18 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.082754048106e+00 1 KSP Residual norm 5.992300533204e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 19 KSP Residual norm 2.890735677419e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.061699839350e+00 ^C[mpiexec@edwards] Sending Ctrl-C to processes as requested and it keeps on going like this. 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. Options used: -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 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 0.000000000000e+00 Linear solve converged due to CONVERGED_ATOL iterations 0 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 2.321914527334e+01 1 KSP Residual norm 5.204444045635e-14 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 9.964395910568e+14 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.573745221888e+00 1 KSP Residual norm 3.612233230950e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 1 KSP Residual norm 2.667315756754e+12 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.053280683677e+00 1 KSP Residual norm 5.770507184021e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 2 KSP Residual norm 3.176176643790e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.016103066894e+00 1 KSP Residual norm 4.558255925705e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 3 KSP Residual norm 1.046851012581e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.095000885520e+00 1 KSP Residual norm 6.356842242710e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 4 KSP Residual norm 4.078940666596e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.244953433655e+00 1 KSP Residual norm 6.745548989971e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 5 KSP Residual norm 1.774320878731e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.176603450320e+00 1 KSP Residual norm 5.385862671655e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 6 KSP Residual norm 8.586271019948e+09 Linear solve converged due to CONVERGED_RTOL iterations 6 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.162170277427e+01 1 KSP Residual norm 2.912041048579e-14 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 9.290138250365e-01 1 KSP Residual norm 2.116243919788e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 3.986822705891e+13 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.573745244651e+00 1 KSP Residual norm 3.715168656229e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 1 KSP Residual norm 1.067205796344e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.053498210605e+00 1 KSP Residual norm 4.526737475809e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 2 KSP Residual norm 1.270173568462e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.017546903473e+00 1 KSP Residual norm 4.593425818353e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 3 KSP Residual norm 4.160360750454e+09 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.097951946288e+00 1 KSP Residual norm 6.461246234077e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 4 KSP Residual norm 1.609903325354e+09 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.250441943430e+00 1 KSP Residual norm 7.157658262436e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 5 KSP Residual norm 6.919994735907e+08 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.169858508137e+00 1 KSP Residual norm 5.046051842799e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 6 KSP Residual norm 3.217042876248e+08 Linear solve converged due to CONVERGED_RTOL iterations 6 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 4.649924240988e-01 1 KSP Residual norm 1.068254961936e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 7.194107998977e-04 1 KSP Residual norm 1.744183339526e-18 Linear solve converged due to CONVERGED_RTOL iterations 1 Residual norms for fieldsplit_1_ solve. 0 KSP Residual norm 4.448127329286e+12 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 7.865130501209e-01 1 KSP Residual norm 5.153218330579e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 1 KSP Residual norm 2.638913823938e+12 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.283369467124e+00 1 KSP Residual norm 7.504488207954e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 2 KSP Residual norm 1.426102610607e+12 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.129662815202e+00 1 KSP Residual norm 5.995004855070e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 3 KSP Residual norm 8.025360882010e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.424860285739e+00 1 KSP Residual norm 8.619123051910e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 4 KSP Residual norm 5.286840519981e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.289665727094e+00 1 KSP Residual norm 7.046395176709e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 5 KSP Residual norm 3.634745308165e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.364708363745e+00 1 KSP Residual norm 7.402799565707e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 6 KSP Residual norm 1.720693722234e+11 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.102971050925e+00 1 KSP Residual norm 6.835965584112e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 7 KSP Residual norm 7.870496175081e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.345118128097e+00 1 KSP Residual norm 7.734345056964e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 8 KSP Residual norm 2.921176674626e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.291212792694e+00 1 KSP Residual norm 7.780332549058e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 9 KSP Residual norm 1.148381817318e+10 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.073046392930e+00 1 KSP Residual norm 6.464513846538e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 10 KSP Residual norm 5.253676669489e+09 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.016140779574e+00 1 KSP Residual norm 5.990033309135e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 11 KSP Residual norm 2.216200526999e+09 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.037961855022e+00 1 KSP Residual norm 6.208761778424e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 12 KSP Residual norm 9.013987504565e+08 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.033050143801e+00 1 KSP Residual norm 5.853852773652e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 13 KSP Residual norm 3.751638125875e+08 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 9.581056747534e-01 1 KSP Residual norm 4.028342105749e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 14 KSP Residual norm 1.372411728307e+08 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 1.117841059532e+00 1 KSP Residual norm 4.848302320534e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 15 KSP Residual norm 5.353181748969e+07 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 8.645767106479e-01 1 KSP Residual norm 4.186764891168e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 16 KSP Residual norm 2.421357309016e+07 Linear solve converged due to CONVERGED_RTOL iterations 16 Residual norms for fieldsplit_0_ solve. 0 KSP Residual norm 6.302684114623e-01 1 KSP Residual norm 2.299313600630e-15 Linear solve converged due to CONVERGED_RTOL iterations 1 Linear solve converged due to CONVERGED_RTOL iterations 2 > Thanks, this looks quite reasonable. I'll try to experiment with it. > > >> > 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. >> >> Okay, there's nothing wrong with that. >> > Thanks for the confirmation. > >
