Great to know there is an example now.
Yes, -pc_fieldsplit_defect_saddle_point worked.
But I would like to understand what I did wrong (and why it didn’t work).
After reading many posted posts, I needed to create PCFieldSplitSetIS to split
the fields.
I set the first N indices to the field phi and the next C indices (starting
from N) to the field lambda.
Please tell me what I did wrong and how I can fix the lines which are now
commented?
PetscErrorCode ierr;
KSP ksp;
PC pc;
KSPCreate(PETSC_COMM_WORLD, &ksp);
KSPSetType(ksp, KSPGMRES);
KSPSetOperators(ksp, A[level], A[level]);
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCFIELDSPLIT);CHKERRQ(ierr);
……..
for (int i = 0; i < N; i++)
vec_zero_field_index[i] = i;
PetscInt * vec_first_field_index = NULL;
PetscMalloc(n * sizeof(PetscInt), &vec_first_field_index);
for (int i = 0; i < C; i++)
vec_first_field_index[i] = N + i;
IS zero_field_isx, first_field_isx;
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, N, vec_zero_field_index,
PETSC_COPY_VALUES, &zero_field_isx));
CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD, C, vec_first_field_index,
PETSC_COPY_VALUES, &first_field_isx));
ierr = PCFieldSplitSetIS(pc,"0",zero_field_isx); CHKERRQ(ierr);
ierr = PCFieldSplitSetIS(pc,"1",first_field_isx); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-ksp_type", "fgmres"); CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-pc_type", "fieldsplit"); CHKERRQ(ierr);
/*
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_ksp_type", "gmres");
CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_0_pc_type", "jacobi");
CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_ksp_type", "preonly");
CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL,"-fieldsplit_1_pc_type", "lu");
CHKERRQ(ierr);*/
ierr = PetscOptionsSetValue(NULL, "-pc_fieldsplit_detect_saddle_point",
NULL);CHKERRQ(ierr);
ierr = PetscOptionsSetValue(NULL, "-options_left", NULL);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
KSPSolve(ksp, b, x);
Best,
Karthik.
From: Matthew Knepley <knep...@gmail.com>
Date: Monday, 24 April 2023 at 17:57
To: Chockalingam, Karthikeyan (STFC,DL,HC) <karthikeyan.chockalin...@stfc.ac.uk>
Cc: petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Mon, Apr 24, 2023 at 10:22 AM Karthikeyan Chockalingam - STFC UKRI
<karthikeyan.chockalin...@stfc.ac.uk<mailto:karthikeyan.chockalin...@stfc.ac.uk>>
wrote:
Hello,
I was able to construct the below K matrix (using submatrices P and P^T), which
is of type MATAIJ
K = [A P^T
P 0]
and solved them using a direct solver.
I modified your example to create either AIJ or Nest matrices, but use the same
assembly code:
https://gitlab.com/petsc/petsc/-/merge_requests/6368
However, I was reading online that this is a saddle point problem and I should
be employing PCFIELDSPLIT.
Since I have one monolithic matrix K, I was not sure how to split the fields.
With this particular matrix, you can use
-pc_fieldsplit_detect_saddle_point
and it will split it automatically.
Thanks,
Matt
Best regards,
Karthik.
From: Chockalingam, Karthikeyan (STFC,DL,HC)
<karthikeyan.chockalin...@stfc.ac.uk<mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Date: Wednesday, 19 April 2023 at 17:52
To: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
I have declared the mapping
ISLocalToGlobalMapping mapping;
ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD, 1, n, nindices,
PETSC_COPY_VALUES, &mapping);
But when I use MatSetValuesLocal(), how do I know the above mapping is employed
because it is not one of the parameters passed to the function?
Thank you.
Kind regards,
Karthik.
From: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>
Date: Tuesday, 18 April 2023 at 16:21
To: Chockalingam, Karthikeyan (STFC,DL,HC)
<karthikeyan.chockalin...@stfc.ac.uk<mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Tue, Apr 18, 2023 at 11:16 AM Karthikeyan Chockalingam - STFC UKRI
<karthikeyan.chockalin...@stfc.ac.uk<mailto:karthikeyan.chockalin...@stfc.ac.uk>>
wrote:
Thank you for your response. I spend some time understanding how
MatSetValuesLocal and ISLocalToGlobalMappingCreate work.
You can look at SNES ex28 where we do this with DMCOMPOSITE.
Q1) Will the matrix K be of type MATMPIAIJ or MATIS?
K = [A P^T
P 0]
I assume MPIAIJ since IS is only used for Neumann-Neumann decompositions.
Q2) Can I use both MatSetValues() to MatSetValuesLocal() to populate K? Since I
have already used MatSetValues() to construct A.
You can, and there would be no changes in serial if K is exactly the upper left
block, but in parallel global indices would change.
Q3) What are the advantages of using MatSetValuesLocal()? Is it that I can
construct P directly using local indies and map the entrees to the global index
in K?
You have a monolithic K, so that you can use sparse direct solvers to check
things. THis is impossible with separate storage.
Q4) I probably don’t have to construct an independent P matrix
You wouldn't in this case.
Thanks,
Matt
Best regards,
Karthik.
From: Matthew Knepley <knep...@gmail.com<mailto:knep...@gmail.com>>
Date: Tuesday, 18 April 2023 at 11:08
To: Chockalingam, Karthikeyan (STFC,DL,HC)
<karthikeyan.chockalin...@stfc.ac.uk<mailto:karthikeyan.chockalin...@stfc.ac.uk>>
Cc: petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>
<petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>>
Subject: Re: [petsc-users] Setting up a matrix for Lagrange multiplier
On Tue, Apr 18, 2023 at 5:24 AM Karthikeyan Chockalingam - STFC UKRI via
petsc-users <petsc-users@mcs.anl.gov<mailto:petsc-users@mcs.anl.gov>> wrote:
Hello,
I'm solving a problem using the Lagrange multiplier, the matrix has the form
K = [A P^T
P 0]
I am familiar with constructing K using MATMPIAIJ. However, I would like to
know if had [A], can I augment it with [P], [P^T] and [0] of type MATMPIAIJ?
Likewise for vectors as well.
Can you please point me to the right resource, if it is a common operation in
PETSc?
You can do this at least 2 ways:
1) Assemble you submatrices directly into the larger matrix by constructing
local-to-global maps for the emplacement. so that you do
not change your assembly code, except to change MatSetValues() to
MatSetValuesLocal(). This is usually preferable.
2) Use MATNEST and VecNEST to put pointers to submatrices and subvectors
directly in.
Thanks,
Matt
Many thanks.
Kind regards,
Karthik.
--
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
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
--
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
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
--
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
https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
