[petsc-users] DMComposite and Matrix construction

[Barry Smith]

> On Jul 31, 2015, at 11:34 AM, Gautam Bisht <gbisht at lbl.gov> wrote:
> 
> I have a followup question on this topic regarding memory usage. 
> 
> For the Jacobian matrix in src/snes/examples/tutorials/ex28.c of problem_type 2,  MAT_NEW_NONZERO_LOCATION_ERR and MAT_NEW_NONZERO_ALLOCATION_ERR are set PETSC_FALSE. If the problem_type 2 is solved multiple times (see the git difference below), would the memory usage keep on increasing? Or, since the non-zero pattern isn't changing over time, there would only be an increase in memory on the first SNESSolve?

  Only the first time through. In fact after your first matrix assembly you could turn those flags back on to generate errors if you start adding new locations the second time (just as a double check that there is not some mistake in your code the generates matrix entries).

  Barry

> 
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> >git diff ex28.c
> diff --git a/src/snes/examples/tutorials/ex28.c b/src/snes/examples/tutorials/ex28.c
> index 98dd6b9..40f94d1 100644
> --- a/src/snes/examples/tutorials/ex28.c
> +++ b/src/snes/examples/tutorials/ex28.c
> @@ -342,6 +342,7 @@ int main(int argc, char *argv[])
>    Mat            B;
>    IS             *isg;
>    PetscBool      view_draw,pass_dm;
> +  PetscInt iter;
>  
>    PetscInitialize(&argc,&argv,0,help);
>    ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-10,1,1,NULL,&dau);CHKERRQ(ierr);
> @@ -434,7 +435,10 @@ int main(int argc, char *argv[])
>         * of splits, but it requires using a DM (perhaps your own implementation). */
>        ierr = SNESSetDM(snes,pack);CHKERRQ(ierr);
>      }
> +    for (iter=1; iter<10000; iter++)  {
>      ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr);
> +    ierr = FormInitial_Coupled(user,X);CHKERRQ(ierr);
> +    }
>      break;
>    }
>    if (view_draw) {ierr = VecView(X,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);}
> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> 
> 
> -Gautam.
> 
> On Wed, Jul 29, 2015 at 11:25 AM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> > On Jul 29, 2015, at 12:58 PM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> >
> > That fixed it, thanks.  Given that I the nonzero pattern will not change throughout the solve, would it still be unreasonable to use this for larger problems?
> 
>    Eventually it will kill you in the first time you assemble the matrix.
> 
>   Barry
> 
> >
> > -gideon
> >
> >> On Jul 29, 2015, at 1:52 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >>
> >>
> >>   To start, immediately after you have created the matrix do what it says call
> >>
> >>> MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)
> >>
> >>  this will allow you to build your matrix and work fine for modest size problems. For larger problems (only when you have your code running and solving the problem you want) you need to add in more preallocation information. But to do it now would be premature optimization.
> >>
> >>   Barry
> >>
> >>
> >>
> >>> On Jul 29, 2015, at 11:51 AM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> >>>
> >>> That’s generating malloc errors:
> >>>
> >>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
> >>> [0]PETSC ERROR: Argument out of range
> >>> [0]PETSC ERROR: New nonzero at (10,0) caused a malloc
> >>> Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check
> >>>
> >>> I suspect this is because DMCreateMatrix is picking a sparsity pattern which is not consistent with what I need.
> >>>
> >>> -gideon
> >>>
> >>>> On Jul 28, 2015, at 10:10 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >>>>
> >>>>
> >>>> DMCreateMatrix() ?
> >>>>
> >>>>
> >>>>> On Jul 28, 2015, at 9:02 PM, Gideon Simpson <gideon.simpson at gmail.com> wrote:
> >>>>>
> >>>>> I’m working with a DMComposite where I have a DMRedundant with 2 parameters, and then a standard DMDACreate with some number of entires that I would like to have distributed.  For concreteness, suppose it is
> >>>>>
> >>>>> DMCompositeCreate(PETSC_COMM_WORLD, &packer);
> >>>>> DMRedundantCreate(PETSC_COMM_WORLD, 0, 2, &p_dm);
> >>>>> DMCompositeAddDM(packer,p_dm);
> >>>>> DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, nx, 1, 1, NULL,&u_dm);
> >>>>> DMCompositeAddDM(packer,u_dm);
> >>>>> DMCreateGlobalVector(packer,&U);
> >>>>>
> >>>>> Now, I would like to construct a matrix for this problem that can be used for computing Jacobians in a nonlinear solve.  Is there a way to get the matrix size to layout in a “useful” way, in the sense that the first process, which owns the two degrees of freedom of p_dm and the first N0 number of the rows of u_dm, controls the corresponding N0+2 rows of the matrix, and analgously for the second process has the next N1 rows of the u_dm vector, and has the next N1 rows of the matrix?
> >>>>>
> >>>>> -gideon
> >>>>>
> >>>>
> >>>
> >>
> >
> 
>