[petsc-users] discontinuous viscosity stokes equation 3D staggered grid
Dave May
dave.mayhem23 at gmail.com
Fri Aug 23 08:47:25 CDT 2013
Yes of course. Just shut all the other sub-solvers off
As Matt says, check check check your GAMG configuration.
-ksp_max_it 1
-pc_fieldsplit_0_ksp_max_it 4
-pc_fieldsplit_1_ksp_max_it 2
-ksp_view
Or if you just want to see GAMG configuration, do this
-ksp_max_it 1
-pc_fieldsplit_0_ksp_max_it 4
-pc_fieldsplit_1_ksp_max_it 2
-pc_fieldsplit_0_ksp_view
For your testing, I'd use FMGRES on A00 as this give you the most
flexibility with the choice of smoother.
-pc_fieldsplit_0_ksp_type fgmres
and then do some seriously heavy smoothing on each level to see if you can
make something converge.
By "heavy" I mean something like
gmres(20)+ILU(0)
On 23 August 2013 15:42, Bishesh Khanal <bisheshkh at gmail.com> wrote:
>
>
>
> On Fri, Aug 23, 2013 at 3:34 PM, Matthew Knepley <knepley at gmail.com>wrote:
>
>> On Fri, Aug 23, 2013 at 8:30 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>
>>>
>>>
>>>
>>> On Fri, Aug 23, 2013 at 3:16 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>
>>>> On Fri, Aug 23, 2013 at 8:01 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>
>>>>>
>>>>>
>>>>>
>>>>> On Fri, Aug 23, 2013 at 2:53 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>
>>>>>> On Fri, Aug 23, 2013 at 7:46 AM, Bishesh Khanal <bisheshkh at gmail.com>wrote:
>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> On Fri, Aug 23, 2013 at 2:33 PM, Matthew Knepley <knepley at gmail.com>wrote:
>>>>>>>
>>>>>>>> On Fri, Aug 23, 2013 at 7:25 AM, Bishesh Khanal <
>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Fri, Aug 23, 2013 at 2:09 PM, Matthew Knepley <
>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>
>>>>>>>>>> On Fri, Aug 23, 2013 at 4:31 AM, Bishesh Khanal <
>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Thanks Matt and Mark for comments in using near null space
>>>>>>>>>>> [question I asked in the thread with subject: *problem
>>>>>>>>>>> (Segmentation voilation) using -pc_type hypre -pc_hypre_type -pilut with
>>>>>>>>>>> multiple nodes in a cluster*].
>>>>>>>>>>> So I understood that I have to set a nearNullSpace to A00 block
>>>>>>>>>>> where the null space correspond to the rigid body motion. I tried it but
>>>>>>>>>>> still the gamg just keeps on iterating and convergence is very very slow. I
>>>>>>>>>>> am not sure what the problem is, right now gamg does not even work for the
>>>>>>>>>>> constant viscosity case.
>>>>>>>>>>> I have set up the following in my code:
>>>>>>>>>>> 1. null space for the whole system A 2. null space for the Schur
>>>>>>>>>>> complement S 3. Near null space for A00 4. a user preconditioner matrix of
>>>>>>>>>>> inverse viscosity in the diagonal for S.
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> If you want to debug solvers, you HAVE to send -ksp_view.
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>> When I use gamg, the -fieldsplit_0_ksp was iterating on and on so
>>>>>>>>> didn't get to the end to get -ksp_view results.
>>>>>>>>> Instead here I have put the -ksp_view output when running the
>>>>>>>>> program with following options: (In this case I get the results)
>>>>>>>>> -pc_type fieldsplit -pc_fieldsplit_type schur
>>>>>>>>> -pc_fieldsplit_dm_splits 0 -pc_fieldsplit_0_fields 0,1,2
>>>>>>>>> -pc_fieldsplit_1_fields 3 -ksp_converged_reason -ksp_view
>>>>>>>>>
>>>>>>>>
>>>>>>>> Okay, that looks fine. Does
>>>>>>>>
>>>>>>>> -fieldsplit_0_pc_type lu
>>>>>>>> - fieldsplit_1_ksp_rtol 1.0e-10
>>>>>>>>
>>>>>>>> converge in one Iterate?
>>>>>>>>
>>>>>>>> What matrix did you attach as the preconditioner matrix for
>>>>>>>> fieldsplit_1_?
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>> I used a diagonal matrix with reciprocal of viscosity values of the
>>>>>>> corresponding cell centers as the preconditioner.
>>>>>>>
>>>>>>> with the options -fieldsplit_0_pc_type lu - fieldsplit_1_ksp_rtol
>>>>>>> 1.0e-10 -fieldsplit_1_ksp_converged_reason -ksp_converged_reason
>>>>>>> I get the following output which means the outer ksp did converge in
>>>>>>> one iterate I guess.
>>>>>>> Linear solve converged due to CONVERGED_RTOL iterations 18
>>>>>>> Linear solve converged due to CONVERGED_RTOL iterations 18
>>>>>>> Linear solve converged due to CONVERGED_RTOL iterations 1
>>>>>>>
>>>>>>
>>>>>> Okay, so A_00 is nonsingular, and the system seems to solve alright.
>>>>>> What do you get for
>>>>>>
>>>>>> -fieldsplit_0_ksp_max_it 30
>>>>>> -fieldsplit_0_pc_type gamg
>>>>>> -fieldsplit_0_ksp_converged_reason
>>>>>> -fieldsplit_1_ksp_converged_reason
>>>>>>
>>>>>>
>>>>>
>>>>> It fieldsplit_0_ does not converge in 30 iterations. It gives:
>>>>> Linear solve converged due to CONVERGED_ATOL iterations 0
>>>>> Linear solve did not converge due to DIVERGED_ITS iterations 30
>>>>>
>>>>> and continues with the same message.
>>>>>
>>>>
>>>> So what would you do? Give up?
>>>>
>>>> No, I don't want to give up :)
>>>
>>>
>>>> -fieldsplit_0_ksp_gmres_restart 200
>>>>
>>>> The idea is to figure out what is going on:
>>>>
>>>> -fieldsplit_0_ksp_monitor_true_residual
>>>>
>>>> I have tried these options before too, the residual is decreasing very
>>> very slowly, but I've not been able to figure out why. (using hypre does
>>> converge although slowly again, but I had problems using hypre with
>>> multiple nodes in a cluster with segmentation fault (we discussed that in
>>> another thread!) )
>>>
>>
>> Put in the Laplacian instead of the operator you have now. It should
>> converge in a few iterates. If not, you have a problem
>> in the specification.
>>
>> If so, put in linear elasticity. If it is slow, you have made a mistake
>> specifiying the near null space. Also, you need to check
>> that the near null space made it to GAMG using the ksp_view output.
>>
>
> Which operator are you referring to ? The one in A00 block ? I'm testing
> currently with the constant viscosity case which means the A00 block has
> \mu div(grad(v)) which is a Laplacian.
> And Is it possible to view the ksp_view output before the solver actually
> converges to check if GAMG took the near null space ?
>
>
>>
>> Matt
>>
>>
>>> e.g a snapshot of the output:
>>>
>>> Residual norms for fieldsplit_0_ solve.
>>> 0 KSP preconditioned resid norm 0.000000000000e+00 true resid norm
>>> 0.000000000000e+00 ||r(i)||/||b|| -nan
>>> Linear solve converged due to CONVERGED_ATOL iterations 0
>>> Residual norms for fieldsplit_0_ solve.
>>> 0 KSP preconditioned resid norm 2.619231455875e-01 true resid norm
>>> 3.637306695895e+02 ||r(i)||/||b|| 1.000000000000e+00
>>> 1 KSP preconditioned resid norm 9.351491725479e-02 true resid norm
>>> 6.013334574957e+01 ||r(i)||/||b|| 1.653238255038e-01
>>> 2 KSP preconditioned resid norm 6.010357491087e-02 true resid norm
>>> 3.664473273769e+01 ||r(i)||/||b|| 1.007468871928e-01
>>> 3 KSP preconditioned resid norm 6.006968012944e-02 true resid norm
>>> 3.696451770148e+01 ||r(i)||/||b|| 1.016260678353e-01
>>> 4 KSP preconditioned resid norm 4.418407037098e-02 true resid norm
>>> 3.184810838034e+01 ||r(i)||/||b|| 8.755959022176e-02
>>> ...
>>> ...
>>> 93 KSP preconditioned resid norm 4.549506047737e-04 true resid norm
>>> 2.877594552685e+00 ||r(i)||/||b|| 7.911333283864e-03
>>> 94 KSP preconditioned resid norm 4.515424416235e-04 true resid norm
>>> 2.875249044668e+00 ||r(i)||/||b|| 7.904884809172e-03
>>> 95 KSP preconditioned resid norm 4.277647876573e-04 true resid norm
>>> 2.830418831358e+00 ||r(i)||/||b|| 7.781633686685e-03
>>> 96 KSP preconditioned resid norm 4.244529173876e-04 true resid norm
>>> 2.807041401408e+00 ||r(i)||/||b|| 7.717362422521e-03
>>> 97 KSP preconditioned resid norm 4.138326570674e-04 true resid norm
>>> 2.793663020386e+00 ||r(i)||/||b|| 7.680581413547e-03
>>> 98 KSP preconditioned resid norm 3.869979433609e-04 true resid norm
>>> 2.715150386650e+00 ||r(i)||/||b|| 7.464727650583e-03
>>> 99 KSP preconditioned resid norm 3.847873979265e-04 true resid norm
>>> 2.706008990336e+00 ||r(i)||/||b|| 7.439595328571e-03
>>>
>>> ....
>>> ....
>>> 294 KSP preconditioned resid norm 1.416482289961e-04 true resid norm
>>> 2.735750748819e+00 ||r(i)||/||b|| 7.521363958412e-03
>>> 295 KSP preconditioned resid norm 1.415389087364e-04 true resid norm
>>> 2.742638608355e+00 ||r(i)||/||b|| 7.540300661064e-03
>>> 296 KSP preconditioned resid norm 1.414967651105e-04 true resid norm
>>> 2.747224243968e+00 ||r(i)||/||b|| 7.552907889424e-03
>>> 297 KSP preconditioned resid norm 1.413843018303e-04 true resid norm
>>> 2.752574248710e+00 ||r(i)||/||b|| 7.567616587891e-03
>>> 298 KSP preconditioned resid norm 1.411747949695e-04 true resid norm
>>> 2.765459647367e+00 ||r(i)||/||b|| 7.603042246859e-03
>>> 299 KSP preconditioned resid norm 1.411609742082e-04 true resid norm
>>> 2.765900464868e+00 ||r(i)||/||b|| 7.604254180683e-03
>>> 300 KSP preconditioned resid norm 1.409844332838e-04 true resid norm
>>> 2.771790506811e+00 ||r(i)||/||b|| 7.620447596402e-03
>>> Linear solve did not converge due to DIVERGED_ITS iterations 300
>>> Residual norms for fieldsplit_0_ solve.
>>> 0 KSP preconditioned resid norm 1.294272083271e-03 true resid norm
>>> 1.776945075651e+00 ||r(i)||/||b|| 1.000000000000e+00
>>> ...
>>> ...
>>>
>>>
>>>
>>>
>>>> Matt
>>>>
>>>>
>>>>>
>>>>>
>>>>>
>>>>>> This is the kind of investigation you msut be comfortable with if you
>>>>>> want to experiment with these solvers.
>>>>>>
>>>>>> Matt
>>>>>>
>>>>>>
>>>>>>>
>>>>>>>>
>>>>>>>> Thanks,
>>>>>>>>
>>>>>>>> Matt
>>>>>>>>
>>>>>>>>
>>>>>>>>> Linear solve converged due to CONVERGED_RTOL iterations 2
>>>>>>>>> KSP Object: 1 MPI processes
>>>>>>>>> type: gmres
>>>>>>>>> GMRES: restart=30, using Classical (unmodified) Gram-Schmidt
>>>>>>>>> Orthogonalization with no iterative refinement
>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50, divergence=10000
>>>>>>>>> left preconditioning
>>>>>>>>> has attached null space
>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>> PC Object: 1 MPI processes
>>>>>>>>> type: fieldsplit
>>>>>>>>> FieldSplit with Schur preconditioner, blocksize = 4,
>>>>>>>>> factorization FULL
>>>>>>>>> Preconditioner for the Schur complement formed from user
>>>>>>>>> provided matrix
>>>>>>>>> Split info:
>>>>>>>>> Split number 0 Fields 0, 1, 2
>>>>>>>>> Split number 1 Fields 3
>>>>>>>>> KSP solver for A00 block
>>>>>>>>> KSP Object: (fieldsplit_0_) 1 MPI processes
>>>>>>>>> type: gmres
>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>> divergence=10000
>>>>>>>>> left preconditioning
>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>> PC Object: (fieldsplit_0_) 1 MPI processes
>>>>>>>>> type: ilu
>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>> 0 levels of fill
>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>> using diagonal shift on blocks to prevent zero pivot
>>>>>>>>> matrix ordering: natural
>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>> Factored matrix follows:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=8232, cols=8232
>>>>>>>>> package used to perform factorization: petsc
>>>>>>>>> total: nonzeros=576000, allocated nonzeros=576000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> using I-node routines: found 2744 nodes, limit
>>>>>>>>> used is 5
>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=8232, cols=8232
>>>>>>>>> total: nonzeros=576000, allocated nonzeros=576000
>>>>>>>>> total number of mallocs used during MatSetValues calls =0
>>>>>>>>> using I-node routines: found 2744 nodes, limit used is
>>>>>>>>> 5
>>>>>>>>> KSP solver for S = A11 - A10 inv(A00) A01
>>>>>>>>> KSP Object: (fieldsplit_1_) 1 MPI processes
>>>>>>>>> type: gmres
>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>> divergence=10000
>>>>>>>>> left preconditioning
>>>>>>>>> has attached null space
>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>> PC Object: (fieldsplit_1_) 1 MPI processes
>>>>>>>>> type: ilu
>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>> 0 levels of fill
>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>> using diagonal shift on blocks to prevent zero pivot
>>>>>>>>> matrix ordering: natural
>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>> Factored matrix follows:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=2744, cols=2744
>>>>>>>>> package used to perform factorization: petsc
>>>>>>>>> total: nonzeros=64000, allocated nonzeros=64000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> not using I-node routines
>>>>>>>>> linear system matrix followed by preconditioner matrix:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: schurcomplement
>>>>>>>>> rows=2744, cols=2744
>>>>>>>>> Schur complement A11 - A10 inv(A00) A01
>>>>>>>>> A11
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=2744, cols=2744
>>>>>>>>> total: nonzeros=64000, allocated nonzeros=64000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> not using I-node routines
>>>>>>>>> A10
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=2744, cols=8232
>>>>>>>>> total: nonzeros=192000, allocated nonzeros=192000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> not using I-node routines
>>>>>>>>> KSP of A00
>>>>>>>>> KSP Object:
>>>>>>>>> (fieldsplit_0_) 1 MPI processes
>>>>>>>>> type: gmres
>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>> divergence=10000
>>>>>>>>> left preconditioning
>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>> PC Object:
>>>>>>>>> (fieldsplit_0_) 1 MPI processes
>>>>>>>>> type: ilu
>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>> 0 levels of fill
>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>> using diagonal shift on blocks to prevent zero
>>>>>>>>> pivot
>>>>>>>>> matrix ordering: natural
>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>> Factored matrix follows:
>>>>>>>>> Matrix Object: 1 MPI
>>>>>>>>> processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=8232, cols=8232
>>>>>>>>> package used to perform factorization:
>>>>>>>>> petsc
>>>>>>>>> total: nonzeros=576000, allocated
>>>>>>>>> nonzeros=576000
>>>>>>>>> total number of mallocs used during
>>>>>>>>> MatSetValues calls =0
>>>>>>>>> using I-node routines: found 2744 nodes,
>>>>>>>>> limit used is 5
>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=8232, cols=8232
>>>>>>>>> total: nonzeros=576000, allocated nonzeros=576000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> using I-node routines: found 2744 nodes, limit
>>>>>>>>> used is 5
>>>>>>>>> A01
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=8232, cols=2744
>>>>>>>>> total: nonzeros=192000, allocated nonzeros=192000
>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>> calls =0
>>>>>>>>> using I-node routines: found 2744 nodes, limit
>>>>>>>>> used is 5
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=2744, cols=2744
>>>>>>>>> total: nonzeros=64000, allocated nonzeros=64000
>>>>>>>>> total number of mallocs used during MatSetValues calls =0
>>>>>>>>> not using I-node routines
>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>> type: seqaij
>>>>>>>>> rows=10976, cols=10976, bs=4
>>>>>>>>> total: nonzeros=1024000, allocated nonzeros=1024000
>>>>>>>>> total number of mallocs used during MatSetValues calls =0
>>>>>>>>> using I-node routines: found 2744 nodes, limit used is 5
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>>> Matt
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> I am testing a small problem with CONSTANT viscosity for grid
>>>>>>>>>>> size of 14^3 with the run time option:
>>>>>>>>>>> -ksp_type gcr -pc_type fieldsplit -pc_fieldsplit_type schur
>>>>>>>>>>> -pc_fieldsplit_dm_splits 0 -pc_fieldsplit_0_fields 0,1,2
>>>>>>>>>>> -pc_fieldsplit_1_fields 3 -ksp_converged_reason -ksp_view
>>>>>>>>>>> -fieldsplit_0_ksp_type gcr -fieldsplit_0_pc_type gamg
>>>>>>>>>>> -fieldsplit_0_ksp_monitor_true_residual -fieldsplit_0_ksp_converged_reason
>>>>>>>>>>> -fieldsplit_1_ksp_monitor_true_residual
>>>>>>>>>>>
>>>>>>>>>>> Here is my relevant code of the solve function:
>>>>>>>>>>> PetscErrorCode ierr;
>>>>>>>>>>> PetscFunctionBeginUser;
>>>>>>>>>>> ierr =
>>>>>>>>>>> DMKSPSetComputeRHS(mDa,computeRHSTaras3D,this);CHKERRQ(ierr);
>>>>>>>>>>> ierr =
>>>>>>>>>>> DMKSPSetComputeOperators(mDa,computeMatrixTaras3D,this);CHKERRQ(ierr);
>>>>>>>>>>> ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr); //mDa
>>>>>>>>>>> with dof = 4, vx,vy,vz and p.
>>>>>>>>>>> ierr =
>>>>>>>>>>> KSPSetNullSpace(mKsp,mNullSpace);CHKERRQ(ierr);//nullSpace for the main
>>>>>>>>>>> system
>>>>>>>>>>> ierr = KSPSetFromOptions(mKsp);CHKERRQ(ierr);
>>>>>>>>>>> ierr = KSPSetUp(mKsp);CHKERRQ(ierr);
>>>>>>>>>>> //register the fieldsplits obtained from options.
>>>>>>>>>>>
>>>>>>>>>>> //Setting up user PC for Schur Complement
>>>>>>>>>>> ierr = KSPGetPC(mKsp,&mPc);CHKERRQ(ierr);
>>>>>>>>>>> ierr =
>>>>>>>>>>> PCFieldSplitSchurPrecondition(mPc,PC_FIELDSPLIT_SCHUR_PRE_USER,mPcForSc);CHKERRQ(ierr);
>>>>>>>>>>>
>>>>>>>>>>> KSP *subKsp;
>>>>>>>>>>> PetscInt subKspPos = 0;
>>>>>>>>>>> //Set up nearNullspace for A00 block.
>>>>>>>>>>> ierr =
>>>>>>>>>>> PCFieldSplitGetSubKSP(mPc,&subKspPos,&subKsp);CHKERRQ(ierr);
>>>>>>>>>>> MatNullSpace rigidBodyModes;
>>>>>>>>>>> Vec coords;
>>>>>>>>>>> ierr = DMGetCoordinates(mDa,&coords);CHKERRQ(ierr);
>>>>>>>>>>> ierr =
>>>>>>>>>>> MatNullSpaceCreateRigidBody(coords,&rigidBodyModes);CHKERRQ(ierr);
>>>>>>>>>>> Mat matA00;
>>>>>>>>>>> ierr =
>>>>>>>>>>> KSPGetOperators(subKsp[0],&matA00,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>> ierr =
>>>>>>>>>>> MatSetNearNullSpace(matA00,rigidBodyModes);CHKERRQ(ierr);
>>>>>>>>>>> ierr = MatNullSpaceDestroy(&rigidBodyModes);CHKERRQ(ierr);
>>>>>>>>>>>
>>>>>>>>>>> //Position 1 => Ksp corresponding to Schur complement S on
>>>>>>>>>>> pressure space
>>>>>>>>>>> subKspPos = 1;
>>>>>>>>>>> ierr =
>>>>>>>>>>> PCFieldSplitGetSubKSP(mPc,&subKspPos,&subKsp);CHKERRQ(ierr);
>>>>>>>>>>> //Set up the null space of constant pressure.
>>>>>>>>>>> ierr = KSPSetNullSpace(subKsp[1],mNullSpaceP);CHKERRQ(ierr);
>>>>>>>>>>> PetscBool isNull;
>>>>>>>>>>> Mat matSc;
>>>>>>>>>>> ierr =
>>>>>>>>>>> KSPGetOperators(subKsp[1],&matSc,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>> ierr = MatNullSpaceTest(mNullSpaceP,matSc,&isNull);
>>>>>>>>>>> if(!isNull)
>>>>>>>>>>> SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"not a valid
>>>>>>>>>>> pressure null space \n");
>>>>>>>>>>> ierr = KSPGetOperators(mKsp,&mA,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>> ierr = MatNullSpaceTest(mNullSpace,mA,&isNull);CHKERRQ(ierr);
>>>>>>>>>>> if(!isNull)
>>>>>>>>>>> SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"not a valid
>>>>>>>>>>> system null space \n");
>>>>>>>>>>>
>>>>>>>>>>> ierr = PetscFree(subKsp);CHKERRQ(ierr);
>>>>>>>>>>> ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>> ierr = KSPGetSolution(mKsp,&mX);CHKERRQ(ierr);
>>>>>>>>>>> ierr = KSPGetRhs(mKsp,&mB);CHKERRQ(ierr);
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> PetscFunctionReturn(0);
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> On Wed, Aug 7, 2013 at 2:15 PM, Matthew Knepley <
>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>
>>>>>>>>>>>> On Wed, Aug 7, 2013 at 7:07 AM, Bishesh Khanal <
>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> On Tue, Aug 6, 2013 at 11:34 PM, Matthew Knepley <
>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>
>>>>>>>>>>>>>> On Tue, Aug 6, 2013 at 10:59 AM, Bishesh Khanal <
>>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> On Tue, Aug 6, 2013 at 4:40 PM, Matthew Knepley <
>>>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> On Tue, Aug 6, 2013 at 8:06 AM, Bishesh Khanal <
>>>>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> On Mon, Aug 5, 2013 at 4:14 PM, Matthew Knepley <
>>>>>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> On Mon, Aug 5, 2013 at 8:48 AM, Bishesh Khanal <
>>>>>>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> On Mon, Aug 5, 2013 at 3:17 PM, Matthew Knepley <
>>>>>>>>>>>>>>>>>>> knepley at gmail.com> wrote:
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> On Mon, Aug 5, 2013 at 7:54 AM, Bishesh Khanal <
>>>>>>>>>>>>>>>>>>>> bisheshkh at gmail.com> wrote:
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>> On Wed, Jul 17, 2013 at 9:48 PM, Jed Brown <
>>>>>>>>>>>>>>>>>>>>> jedbrown at mcs.anl.gov> wrote:
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>> Bishesh Khanal <bisheshkh at gmail.com> writes:
>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>> > Now, I implemented two different approaches, each
>>>>>>>>>>>>>>>>>>>>>> for both 2D and 3D, in
>>>>>>>>>>>>>>>>>>>>>> > MATLAB. It works for the smaller sizes but I have
>>>>>>>>>>>>>>>>>>>>>> problems solving it for
>>>>>>>>>>>>>>>>>>>>>> > the problem size I need (250^3 grid size).
>>>>>>>>>>>>>>>>>>>>>> > I use staggered grid with p on cell centers, and
>>>>>>>>>>>>>>>>>>>>>> components of v on cell
>>>>>>>>>>>>>>>>>>>>>> > faces. Similar split up of K to cell center and
>>>>>>>>>>>>>>>>>>>>>> faces to account for the
>>>>>>>>>>>>>>>>>>>>>> > variable viscosity case)
>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>> Okay, you're using a staggered-grid finite difference
>>>>>>>>>>>>>>>>>>>>>> discretization of
>>>>>>>>>>>>>>>>>>>>>> variable-viscosity Stokes. This is a common problem
>>>>>>>>>>>>>>>>>>>>>> and I recommend
>>>>>>>>>>>>>>>>>>>>>> starting with PCFieldSplit with Schur complement
>>>>>>>>>>>>>>>>>>>>>> reduction (make that
>>>>>>>>>>>>>>>>>>>>>> work first, then switch to block preconditioner).
>>>>>>>>>>>>>>>>>>>>>> You can use PCLSC or
>>>>>>>>>>>>>>>>>>>>>> (probably better for you), assemble a preconditioning
>>>>>>>>>>>>>>>>>>>>>> matrix containing
>>>>>>>>>>>>>>>>>>>>>> the inverse viscosity in the pressure-pressure block.
>>>>>>>>>>>>>>>>>>>>>> This diagonal
>>>>>>>>>>>>>>>>>>>>>> matrix is a spectrally equivalent (or nearly so,
>>>>>>>>>>>>>>>>>>>>>> depending on
>>>>>>>>>>>>>>>>>>>>>> discretization) approximation of the Schur
>>>>>>>>>>>>>>>>>>>>>> complement. The velocity
>>>>>>>>>>>>>>>>>>>>>> block can be solved with algebraic multigrid. Read
>>>>>>>>>>>>>>>>>>>>>> the PCFieldSplit
>>>>>>>>>>>>>>>>>>>>>> docs (follow papers as appropriate) and let us know
>>>>>>>>>>>>>>>>>>>>>> if you get stuck.
>>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>> I was trying to assemble the inverse viscosity
>>>>>>>>>>>>>>>>>>>>> diagonal matrix to use as the preconditioner for the Schur complement solve
>>>>>>>>>>>>>>>>>>>>> step as you suggested. I've few questions about the ways to implement this
>>>>>>>>>>>>>>>>>>>>> in Petsc:
>>>>>>>>>>>>>>>>>>>>> A naive approach that I can think of would be to
>>>>>>>>>>>>>>>>>>>>> create a vector with its components as reciprocal viscosities of the cell
>>>>>>>>>>>>>>>>>>>>> centers corresponding to the pressure variables, and then create a diagonal
>>>>>>>>>>>>>>>>>>>>> matrix from this vector. However I'm not sure about:
>>>>>>>>>>>>>>>>>>>>> How can I make this matrix, (say S_p) compatible to
>>>>>>>>>>>>>>>>>>>>> the Petsc distribution of the different rows of the main system matrix over
>>>>>>>>>>>>>>>>>>>>> different processors ? The main matrix was created using the DMDA structure
>>>>>>>>>>>>>>>>>>>>> with 4 dof as explained before.
>>>>>>>>>>>>>>>>>>>>> The main matrix correspond to the DMDA with 4 dofs but
>>>>>>>>>>>>>>>>>>>>> for the S_p matrix would correspond to only pressure space. Should the
>>>>>>>>>>>>>>>>>>>>> distribution of the rows of S_p among different processor not correspond to
>>>>>>>>>>>>>>>>>>>>> the distribution of the rhs vector, say h' if it is solving for p with Sp =
>>>>>>>>>>>>>>>>>>>>> h' where S = A11 inv(A00) A01 ?
>>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> PETSc distributed vertices, not dofs, so it never
>>>>>>>>>>>>>>>>>>>> breaks blocks. The P distribution is the same as the entire problem divided
>>>>>>>>>>>>>>>>>>>> by 4.
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> Thanks Matt. So if I create a new DMDA with same grid
>>>>>>>>>>>>>>>>>>> size but with dof=1 instead of 4, the vertices for this new DMDA will be
>>>>>>>>>>>>>>>>>>> identically distributed as for the original DMDA ? Or should I inform PETSc
>>>>>>>>>>>>>>>>>>> by calling a particular function to make these two DMDA have identical
>>>>>>>>>>>>>>>>>>> distribution of the vertices ?
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> Yes.
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>> Even then I think there might be a problem due to the
>>>>>>>>>>>>>>>>>>> presence of "fictitious pressure vertices". The system matrix (A) contains
>>>>>>>>>>>>>>>>>>> an identity corresponding to these fictitious pressure nodes, thus when
>>>>>>>>>>>>>>>>>>> using a -pc_fieldsplit_detect_saddle_point, will detect a A11 zero block of
>>>>>>>>>>>>>>>>>>> size that correspond to only non-fictitious P-nodes. So the preconditioner
>>>>>>>>>>>>>>>>>>> S_p for the Schur complement outer solve with Sp = h' will also need to
>>>>>>>>>>>>>>>>>>> correspond to only the non-fictitious P-nodes. This means its size does not
>>>>>>>>>>>>>>>>>>> directly correspond to the DMDA grid defined for the original problem.
>>>>>>>>>>>>>>>>>>> Could you please suggest an efficient way of assembling this S_p matrix ?
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> Don't use detect_saddle, but split it by fields
>>>>>>>>>>>>>>>>>> -pc_fieldsplit_0_fields 0,1,2 -pc_fieldsplit_1_fields 4
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>> How can I set this split in the code itself without giving
>>>>>>>>>>>>>>>>> it as a command line option when the system matrix is assembled from the
>>>>>>>>>>>>>>>>> DMDA for the whole system with 4 dofs. (i.e. *without*using the DMComposite or
>>>>>>>>>>>>>>>>> *without* using the nested block matrices to assemble
>>>>>>>>>>>>>>>>> different blocks separately and then combine them together).
>>>>>>>>>>>>>>>>> I need the split to get access to the fieldsplit_1_ksp in
>>>>>>>>>>>>>>>>> my code, because not using detect_saddle_point means I cannot use
>>>>>>>>>>>>>>>>> -fieldsplit_1_ksp_constant_null_space due to the presence of identity for
>>>>>>>>>>>>>>>>> the fictitious pressure nodes present in the fieldsplit_1_ block. I need to
>>>>>>>>>>>>>>>>> use PCFieldSplitGetSubKsp() so that I can set proper null-space basis.
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> This is currently a real problem with the DMDA. In the
>>>>>>>>>>>>>>>> unstructured case, where we always need specialized spaces, you can
>>>>>>>>>>>>>>>> use something like
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> PetscObject pressure;
>>>>>>>>>>>>>>>> MatNullSpace nullSpacePres;
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> ierr = DMGetField(dm, 1, &pressure);CHKERRQ(ierr);
>>>>>>>>>>>>>>>> ierr = MatNullSpaceCreate(PetscObjectComm(pressure),
>>>>>>>>>>>>>>>> PETSC_TRUE, 0, NULL, &nullSpacePres);CHKERRQ(ierr);
>>>>>>>>>>>>>>>> ierr = PetscObjectCompose(pressure, "nullspace",
>>>>>>>>>>>>>>>> (PetscObject) nullSpacePres);CHKERRQ(ierr);
>>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>>> MatNullSpaceDestroy(&nullSpacePres);CHKERRQ(ierr);
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> and then DMGetSubDM() uses this information to attach the
>>>>>>>>>>>>>>>> null space to the IS that is created using the information in the
>>>>>>>>>>>>>>>> PetscSection.
>>>>>>>>>>>>>>>> If you use a PetscSection to set the data layout over the
>>>>>>>>>>>>>>>> DMDA, I think this works correctly, but this has not been tested at all and
>>>>>>>>>>>>>>>> is very
>>>>>>>>>>>>>>>> new code. Eventually, I think we want all DMs to use this
>>>>>>>>>>>>>>>> mechanism, but we are still working it out.
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> Currently I do not use PetscSection. If this makes a cleaner
>>>>>>>>>>>>>>> approach, I'd try it too but may a bit later (right now I'd like test my
>>>>>>>>>>>>>>> model with a quickfix even if it means a little dirty code!)
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Bottom line: For custom null spaces using the default
>>>>>>>>>>>>>>>> layout in DMDA, you need to take apart the PCFIELDSPLIT after it has been
>>>>>>>>>>>>>>>> setup,
>>>>>>>>>>>>>>>> which is somewhat subtle. You need to call KSPSetUp() and
>>>>>>>>>>>>>>>> then reach in and get the PC, and the subKSPs. I don't like this at all,
>>>>>>>>>>>>>>>> but we
>>>>>>>>>>>>>>>> have not reorganized that code (which could be very simple
>>>>>>>>>>>>>>>> and inflexible since its very structured).
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> So I tried to get this approach working but I could not
>>>>>>>>>>>>>>> succeed and encountered some errors. Here is a code snippet:
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> //mDa is the DMDA that describes the whole grid with all 4
>>>>>>>>>>>>>>> dofs (3 velocity components and 1 pressure comp.)
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> DMKSPSetComputeRHS(mDa,computeRHSTaras3D,this);CHKERRQ(ierr);
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> DMKSPSetComputeOperators(mDa,computeMatrixTaras3D,this);CHKERRQ(ierr);
>>>>>>>>>>>>>>> ierr = KSPSetDM(mKsp,mDa);CHKERRQ(ierr);
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> KSPSetNullSpace(mKsp,mNullSpaceSystem);CHKERRQ(ierr); //I've the
>>>>>>>>>>>>>>> mNullSpaceSystem based on mDa, that contains a null space basis for the
>>>>>>>>>>>>>>> complete system.
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> KSPSetFromOptions(mKsp);CHKERRQ(ierr);
>>>>>>>>>>>>>>> //This I expect would register these options I give:-pc_type fieldsplit
>>>>>>>>>>>>>>> -pc_fieldsplit_type schur -pc_fieldsplit_0_fields 0,1,2
>>>>>>>>>>>>>>> //-pc_fieldsplit_1_fields 3
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ierr = KSPSetUp(mKsp);CHKERRQ(ierr);
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ierr = KSPGetPC(mKsp,&mPcOuter); //Now get the PC
>>>>>>>>>>>>>>> that was obtained from the options (fieldsplit)
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> PCFieldSplitSchurPrecondition(mPcOuter,PC_FIELDSPLIT_SCHUR_PRE_USER,mPcForSc);CHKERRQ(ierr);
>>>>>>>>>>>>>>> //I have created the matrix mPcForSc using a DMDA with identical //size to
>>>>>>>>>>>>>>> mDa but with dof=1 corresponding to the pressure nodes (say mDaPressure).
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ierr = PCSetUp(mPcOuter);CHKERRQ(ierr);
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> KSP *kspSchur;
>>>>>>>>>>>>>>> PetscInt kspSchurPos = 1;
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> PCFieldSplitGetSubKSP(mPcOuter,&kspSchurPos,&kspSchur);CHKERRQ(ierr);
>>>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>>>> KSPSetNullSpace(kspSchur[1],mNullSpacePressure);CHKERRQ(ierr);
>>>>>>>>>>>>>>> //The null space is the one that correspond to only pressure nodes, created
>>>>>>>>>>>>>>> using the mDaPressure.
>>>>>>>>>>>>>>> ierr = PetscFree(kspSchur);CHKERRQ(ierr);
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Sorry, you need to return to the old DMDA behavior, so you
>>>>>>>>>>>>>> want
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> -pc_fieldsplit_dm_splits 0
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> Thanks, with this it seems I can attach the null space
>>>>>>>>>>>>> properly, but I have a question regarding whether the Schur complement ksp
>>>>>>>>>>>>> solver is actually using the preconditioner matrix I provide.
>>>>>>>>>>>>> When using -ksp_view, the outer level pc object of type
>>>>>>>>>>>>> fieldsplit does report that: "Preconditioner for the Schur complement
>>>>>>>>>>>>> formed from user provided matrix", but in the KSP solver for Schur
>>>>>>>>>>>>> complement S, the pc object (fieldsplit_1_) is of type ilu and doesn't say
>>>>>>>>>>>>> that it is using the matrix I provide. Am I missing something here ?
>>>>>>>>>>>>> Below are the relevant commented code snippet and the output
>>>>>>>>>>>>> of the -ksp_view
>>>>>>>>>>>>> (The options I 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 -ksp_converged_reason -ksp_view )
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> If ILU does not error, it means it is using your matrix,
>>>>>>>>>>>> because the Schur complement matrix cannot be factored, and FS says it is
>>>>>>>>>>>> using your matrix.
>>>>>>>>>>>>
>>>>>>>>>>>> Matt
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>> Code snippet:
>>>>>>>>>>>>> ierr = KSPSetNullSpace(mKsp,mNullSpaceSystem);CHKERRQ(ierr);
>>>>>>>>>>>>> //The nullspace for the whole system
>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>> KSPSetFromOptions(mKsp);CHKERRQ(ierr);
>>>>>>>>>>>>> ierr = KSPSetUp(mKsp);CHKERRQ(ierr);
>>>>>>>>>>>>> //Set up mKsp with the options provided with fieldsplit and the fields
>>>>>>>>>>>>> associated with the two splits.
>>>>>>>>>>>>>
>>>>>>>>>>>>> ierr = KSPGetPC(mKsp,&mPcOuter);CHKERRQ(ierr);
>>>>>>>>>>>>> //Get the fieldsplit pc set up from the options
>>>>>>>>>>>>>
>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>> PCFieldSplitSchurPrecondition(mPcOuter,PC_FIELDSPLIT_SCHUR_PRE_USER,mPcForSc);CHKERRQ(ierr);
>>>>>>>>>>>>> //Use mPcForSc as the preconditioner for Schur Complement
>>>>>>>>>>>>>
>>>>>>>>>>>>> KSP *kspSchur;
>>>>>>>>>>>>> PetscInt kspSchurPos = 1;
>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>> PCFieldSplitGetSubKSP(mPcOuter,&kspSchurPos,&kspSchur);CHKERRQ(ierr);
>>>>>>>>>>>>> ierr =
>>>>>>>>>>>>> KSPSetNullSpace(kspSchur[1],mNullSpacePressure);CHKERRQ(ierr);
>>>>>>>>>>>>> //Attach the null-space for the Schur complement ksp solver.
>>>>>>>>>>>>> ierr = PetscFree(kspSchur);CHKERRQ(ierr);
>>>>>>>>>>>>>
>>>>>>>>>>>>> ierr = KSPSolve(mKsp,NULL,NULL);CHKERRQ(ierr);
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>> the output of the -ksp_view
>>>>>>>>>>>>> KSP Object: 1 MPI processes
>>>>>>>>>>>>> type: gmres
>>>>>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50, divergence=10000
>>>>>>>>>>>>> left preconditioning
>>>>>>>>>>>>> has attached null space
>>>>>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>>>>>> PC Object: 1 MPI processes
>>>>>>>>>>>>> type: fieldsplit
>>>>>>>>>>>>> FieldSplit with Schur preconditioner, blocksize = 4,
>>>>>>>>>>>>> factorization FULL
>>>>>>>>>>>>> Preconditioner for the Schur complement formed from user
>>>>>>>>>>>>> provided matrix
>>>>>>>>>>>>> Split info:
>>>>>>>>>>>>> Split number 0 Fields 0, 1, 2
>>>>>>>>>>>>> Split number 1 Fields 3
>>>>>>>>>>>>> KSP solver for A00 block
>>>>>>>>>>>>> KSP Object: (fieldsplit_0_) 1 MPI processes
>>>>>>>>>>>>> type: gmres
>>>>>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>>>>>> divergence=10000
>>>>>>>>>>>>> left preconditioning
>>>>>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>>>>>> PC Object: (fieldsplit_0_) 1 MPI processes
>>>>>>>>>>>>> type: ilu
>>>>>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>>>>>> 0 levels of fill
>>>>>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>>>>>> using diagonal shift on blocks to prevent zero pivot
>>>>>>>>>>>>> matrix ordering: natural
>>>>>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>>>>>> Factored matrix follows:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2187, cols=2187
>>>>>>>>>>>>> package used to perform factorization: petsc
>>>>>>>>>>>>> total: nonzeros=140625, allocated
>>>>>>>>>>>>> nonzeros=140625
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> using I-node routines: found 729 nodes,
>>>>>>>>>>>>> limit used is 5
>>>>>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2187, cols=2187
>>>>>>>>>>>>> total: nonzeros=140625, allocated nonzeros=140625
>>>>>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>>>>>> calls =0
>>>>>>>>>>>>> using I-node routines: found 729 nodes, limit used
>>>>>>>>>>>>> is 5
>>>>>>>>>>>>> KSP solver for S = A11 - A10 inv(A00) A01
>>>>>>>>>>>>> KSP Object: (fieldsplit_1_) 1 MPI processes
>>>>>>>>>>>>> type: gmres
>>>>>>>>>>>>> GMRES: restart=30, using Classical (unmodified)
>>>>>>>>>>>>> Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>>>>>> divergence=10000
>>>>>>>>>>>>> left preconditioning
>>>>>>>>>>>>> has attached null space
>>>>>>>>>>>>> using PRECONDITIONED norm type for convergence test
>>>>>>>>>>>>> PC Object: (fieldsplit_1_) 1 MPI processes
>>>>>>>>>>>>> type: ilu
>>>>>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>>>>>> 0 levels of fill
>>>>>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>>>>>> using diagonal shift on blocks to prevent zero pivot
>>>>>>>>>>>>> matrix ordering: natural
>>>>>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>>>>>> Factored matrix follows:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=729, cols=729
>>>>>>>>>>>>> package used to perform factorization: petsc
>>>>>>>>>>>>> total: nonzeros=15625, allocated nonzeros=15625
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> not using I-node routines
>>>>>>>>>>>>> linear system matrix followed by preconditioner matrix:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: schurcomplement
>>>>>>>>>>>>> rows=729, cols=729
>>>>>>>>>>>>> Schur complement A11 - A10 inv(A00) A01
>>>>>>>>>>>>> A11
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=729, cols=729
>>>>>>>>>>>>> total: nonzeros=15625, allocated nonzeros=15625
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> not using I-node routines
>>>>>>>>>>>>> A10
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=729, cols=2187
>>>>>>>>>>>>> total: nonzeros=46875, allocated nonzeros=46875
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> not using I-node routines
>>>>>>>>>>>>> KSP of A00
>>>>>>>>>>>>> KSP Object:
>>>>>>>>>>>>> (fieldsplit_0_) 1 MPI processes
>>>>>>>>>>>>> type: gmres
>>>>>>>>>>>>> GMRES: restart=30, using Classical
>>>>>>>>>>>>> (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
>>>>>>>>>>>>> GMRES: happy breakdown tolerance 1e-30
>>>>>>>>>>>>> maximum iterations=10000, initial guess is zero
>>>>>>>>>>>>> tolerances: relative=1e-05, absolute=1e-50,
>>>>>>>>>>>>> divergence=10000
>>>>>>>>>>>>> left preconditioning
>>>>>>>>>>>>> using PRECONDITIONED norm type for convergence
>>>>>>>>>>>>> test
>>>>>>>>>>>>> PC Object:
>>>>>>>>>>>>> (fieldsplit_0_) 1 MPI processes
>>>>>>>>>>>>> type: ilu
>>>>>>>>>>>>> ILU: out-of-place factorization
>>>>>>>>>>>>> 0 levels of fill
>>>>>>>>>>>>> tolerance for zero pivot 2.22045e-14
>>>>>>>>>>>>> using diagonal shift on blocks to prevent
>>>>>>>>>>>>> zero pivot
>>>>>>>>>>>>> matrix ordering: natural
>>>>>>>>>>>>> factor fill ratio given 1, needed 1
>>>>>>>>>>>>> Factored matrix follows:
>>>>>>>>>>>>> Matrix Object: 1
>>>>>>>>>>>>> MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2187, cols=2187
>>>>>>>>>>>>> package used to perform factorization:
>>>>>>>>>>>>> petsc
>>>>>>>>>>>>> total: nonzeros=140625, allocated
>>>>>>>>>>>>> nonzeros=140625
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> using I-node routines: found 729
>>>>>>>>>>>>> nodes, limit used is 5
>>>>>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2187, cols=2187
>>>>>>>>>>>>> total: nonzeros=140625, allocated
>>>>>>>>>>>>> nonzeros=140625
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> using I-node routines: found 729 nodes,
>>>>>>>>>>>>> limit used is 5
>>>>>>>>>>>>> A01
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2187, cols=729
>>>>>>>>>>>>> total: nonzeros=46875, allocated nonzeros=46875
>>>>>>>>>>>>> total number of mallocs used during
>>>>>>>>>>>>> MatSetValues calls =0
>>>>>>>>>>>>> using I-node routines: found 729 nodes,
>>>>>>>>>>>>> limit used is 5
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=729, cols=729
>>>>>>>>>>>>> total: nonzeros=15625, allocated nonzeros=15625
>>>>>>>>>>>>> total number of mallocs used during MatSetValues
>>>>>>>>>>>>> calls =0
>>>>>>>>>>>>> not using I-node routines
>>>>>>>>>>>>> linear system matrix = precond matrix:
>>>>>>>>>>>>> Matrix Object: 1 MPI processes
>>>>>>>>>>>>> type: seqaij
>>>>>>>>>>>>> rows=2916, cols=2916, bs=4
>>>>>>>>>>>>> total: nonzeros=250000, allocated nonzeros=250000
>>>>>>>>>>>>> total number of mallocs used during MatSetValues calls =0
>>>>>>>>>>>>> using I-node routines: found 729 nodes, limit used is 5
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> or
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> PCFieldSplitSetDMSplits(pc, PETSC_FALSE)
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Thanks,
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>> The errors I get when running with options: -pc_type
>>>>>>>>>>>>>>> fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_0_fields 0,1,2
>>>>>>>>>>>>>>> -pc_fieldsplit_1_fields 3
>>>>>>>>>>>>>>> [0]PETSC ERROR: --------------------- Error Message
>>>>>>>>>>>>>>> ------------------------------------
>>>>>>>>>>>>>>> [0]PETSC ERROR: No support for this operation for this
>>>>>>>>>>>>>>> object type!
>>>>>>>>>>>>>>> [0]PETSC ERROR: Support only implemented for 2d!
>>>>>>>>>>>>>>> [0]PETSC ERROR:
>>>>>>>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>>>>>>>> [0]PETSC ERROR: Petsc Release Version 3.4.2, Jul, 02, 2013
>>>>>>>>>>>>>>> [0]PETSC ERROR: See docs/changes/index.html for recent
>>>>>>>>>>>>>>> updates.
>>>>>>>>>>>>>>> [0]PETSC ERROR: See docs/faq.html for hints about trouble
>>>>>>>>>>>>>>> shooting.
>>>>>>>>>>>>>>> [0]PETSC ERROR: See docs/index.html for manual pages.
>>>>>>>>>>>>>>> [0]PETSC ERROR:
>>>>>>>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>>>>>>>> [0]PETSC ERROR: src/AdLemMain on a arch-linux2-cxx-debug
>>>>>>>>>>>>>>> named edwards by bkhanal Tue Aug 6 17:35:30 2013
>>>>>>>>>>>>>>> [0]PETSC ERROR: Libraries linked from
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/arch-linux2-cxx-debug/lib
>>>>>>>>>>>>>>> [0]PETSC ERROR: Configure run at Fri Jul 19 14:25:01 2013
>>>>>>>>>>>>>>> [0]PETSC ERROR: Configure options --with-cc=gcc
>>>>>>>>>>>>>>> --with-fc=g77 --with-cxx=g++ --download-f-blas-lapack=1 --download-mpich=1
>>>>>>>>>>>>>>> -with-clanguage=cxx --download-hypre=1
>>>>>>>>>>>>>>> [0]PETSC ERROR:
>>>>>>>>>>>>>>> ------------------------------------------------------------------------
>>>>>>>>>>>>>>> [0]PETSC ERROR: DMCreateSubDM_DA() line 188 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/dm/impls/da/dacreate.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: DMCreateSubDM() line 1267 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/dm/interface/dm.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: PCFieldSplitSetDefaults() line 337 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: PCSetUp_FieldSplit() line 458 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/impls/fieldsplit/fieldsplit.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: PCSetUp() line 890 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/pc/interface/precon.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: KSPSetUp() line 278 in
>>>>>>>>>>>>>>> /home/bkhanal/Documents/softwares/petsc-3.4.2/src/ksp/ksp/interface/itfunc.c
>>>>>>>>>>>>>>> [0]PETSC ERROR: solveModel() line 181 in
>>>>>>>>>>>>>>> "unknowndirectory/"/user/bkhanal/home/works/AdLemModel/src/PetscAdLemTaras3D.cxx
>>>>>>>>>>>>>>> WARNING! There are options you set that were not used!
>>>>>>>>>>>>>>> WARNING! could be spelling mistake, etc!
>>>>>>>>>>>>>>> Option left: name:-pc_fieldsplit_1_fields value: 3
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> Matt
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>> --
>>>>>>>>>>>>>>>>>>>> What most experimenters take for granted before they
>>>>>>>>>>>>>>>>>>>> begin their experiments is infinitely more interesting than any results to
>>>>>>>>>>>>>>>>>>>> which their experiments lead.
>>>>>>>>>>>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>> --
>>>>>>>>>>>>>>>>>> What most experimenters take for granted before they
>>>>>>>>>>>>>>>>>> begin their experiments is infinitely more interesting than any results to
>>>>>>>>>>>>>>>>>> which their experiments lead.
>>>>>>>>>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>> --
>>>>>>>>>>>>>>>> What most experimenters take for granted before they begin
>>>>>>>>>>>>>>>> their experiments is infinitely more interesting than any results to which
>>>>>>>>>>>>>>>> their experiments lead.
>>>>>>>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>>
>>>>>>>>>>>>>> --
>>>>>>>>>>>>>> What most experimenters take for granted before they begin
>>>>>>>>>>>>>> their experiments is infinitely more interesting than any results to which
>>>>>>>>>>>>>> their experiments lead.
>>>>>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>>
>>>>>>>>>>>> --
>>>>>>>>>>>> What most experimenters take for granted before they begin
>>>>>>>>>>>> their experiments is infinitely more interesting than any results to which
>>>>>>>>>>>> their experiments lead.
>>>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>>>> experiments lead.
>>>>>>>>>> -- Norbert Wiener
>>>>>>>>>>
>>>>>>>>>
>>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> --
>>>>>>>> What most experimenters take for granted before they begin their
>>>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>>>> experiments lead.
>>>>>>>> -- Norbert Wiener
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>
>>>>>>
>>>>>> --
>>>>>> What most experimenters take for granted before they begin their
>>>>>> experiments is infinitely more interesting than any results to which their
>>>>>> experiments lead.
>>>>>> -- Norbert Wiener
>>>>>>
>>>>>
>>>>>
>>>>
>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which their
>>>> experiments lead.
>>>> -- Norbert Wiener
>>>>
>>>
>>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their
>> experiments is infinitely more interesting than any results to which their
>> experiments lead.
>> -- Norbert Wiener
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130823/e95b296e/attachment-0001.html>
More information about the petsc-users
mailing list