[petsc-users] solving stokes-like equation in 3D staggered grid in irregular domain with petscsection
Bishesh Khanal
bisheshkh at gmail.com
Tue Dec 3 08:29:25 CST 2013
On Thu, Nov 14, 2013 at 1:05 PM, Bishesh Khanal <bisheshkh at gmail.com> wrote:
>
>
>
> On Wed, Nov 13, 2013 at 10:06 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
>
>> Bishesh Khanal <bisheshkh at gmail.com> 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 at 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.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131203/c5d0ae05/attachment-0001.html>
More information about the petsc-users
mailing list