On Mon, Apr 24, 2023 at 1:37 PM Karthikeyan Chockalingam - STFC UKRI < karthikeyan.chockalin...@stfc.ac.uk> wrote:
> 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). > Can you show the error? Or give something self-contained that I can run. > After reading many posted posts, I needed to create PCFieldSplitSetIS to > split the fields. > Yes. I would give NULL for the name here. > I set the first N indices to the field phi and the next C indices > (starting from N) to the field lambda. > This looks fine. However, these are global indices, so this would not work in parallel. You would have to offset by the first index on this process. Thanks, Matt > 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> 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> > *Date: *Wednesday, 19 April 2023 at 17:52 > *To: *Matthew Knepley <knep...@gmail.com> > *Cc: *petsc-users@mcs.anl.gov <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> > *Date: *Tuesday, 18 April 2023 at 16:21 > *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 Tue, Apr 18, 2023 at 11:16 AM Karthikeyan Chockalingam - STFC UKRI < > 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> > *Date: *Tuesday, 18 April 2023 at 11:08 > *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 Tue, Apr 18, 2023 at 5:24 AM Karthikeyan Chockalingam - STFC UKRI via > petsc-users <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/> > -- 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/>