On Wed, Mar 31, 2021 at 10:15 AM Francesco Brarda <[email protected]> wrote:

Thank you for your advices. I wrote what seems to me a very basic code, but I got this error when I run it with more than 1 processor: Clearly the result 299. is wrong but I do not understand what am doing wrong. With 1 processor it works fine.

My guess is that you do VecGetArray() and index the array using global indices rather than local indices, because there memory corruption with a Vec array. Thanks, Matt

steps 150, ftime 15. Vec Object: 2 MPI processes type: mpi Process [0] 16.5613 2.91405 Process [1] 299. [0]PETSC ERROR: PetscTrFreeDefault() called from VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [0]PETSC ERROR: Block [id=0(16)] at address 0x15812a0 is corrupted (probably write past end of array) [0]PETSC ERROR: Block allocated in VecCreate_MPI_Private() line 514 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pbvec.c [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Memory corruption: https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind [0]PETSC ERROR: Corrupted memory [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Wed Mar 31 16:05:22 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 PetscTrFreeDefault() line 310 in /home/fbrarda/petsc/src/sys/memory/mtr.c [0]PETSC ERROR: #2 VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [0]PETSC ERROR: #3 VecDestroy() line 396 in /home/fbrarda/petsc/src/vec/vec/interface/vector.c [0]PETSC ERROR: #4 SNESLineSearchReset() line 284 in /home/fbrarda/petsc/src/snes/linesearch/interface/linesearch.c [0]PETSC ERROR: #5 SNESReset() line 3229 in /home/fbrarda/petsc/src/snes/interface/snes.c [0]PETSC ERROR: #6 TSReset() line 2800 in /home/fbrarda/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #7 TSDestroy() line 2856 in /home/fbrarda/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #8 main() line 256 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] application called MPI_Abort(MPI_COMM_SELF, 256001) - process 0 [1]PETSC ERROR: PetscTrFreeDefault() called from VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [1]PETSC ERROR: Block [id=0(16)] at address 0xbd9520 is corrupted (probably write past end of array) [1]PETSC ERROR: Block allocated in VecCreate_MPI_Private() line 514 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pbvec.c [1]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [1]PETSC ERROR: Memory corruption: https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind [1]PETSC ERROR: Corrupted memory [1]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [1]PETSC ERROR: Petsc Release Version 3.14.4, unknown [1]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Wed Mar 31 16:05:22 2021 [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [1]PETSC ERROR: #1 PetscTrFreeDefault() line 310 in /home/fbrarda/petsc/src/sys/memory/mtr.c [1]PETSC ERROR: #2 VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [1]PETSC ERROR: #3 VecDestroy() line 396 in /home/fbrarda/petsc/src/vec/vec/interface/vector.c [1]PETSC ERROR: #4 TSReset() line 2806 in /home/fbrarda/petsc/src/ts/interface/ts.c [1]PETSC ERROR: #5 TSDestroy() line 2856 in /home/fbrarda/petsc/src/ts/interface/ts.c [1]PETSC ERROR: #6 main() line 256 in par_sir_model.c [1]PETSC ERROR: No PETSc Option Table entries [1]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] application called MPI_Abort(MPI_COMM_SELF, 256001) - process 0

Il giorno 31 mar 2021, alle ore 12:06, Matthew Knepley <[email protected]> ha scritto:

On Wed, Mar 31, 2021 at 3:54 AM Francesco Brarda < [email protected]> wrote: Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

You do not need one. We use it in examples because it makes it easy to create the data.

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Since you are producing all the Jacobian values, it is whatever is easier in your code I think.

THanks,

Matt

Best, Francesco

-- 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/

-- 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/>

Right now this is only a toy example. I do not expect to see any actual improvement over the number of processors, or should I?

Il giorno 31 mar 2021, alle ore 16:43, Stefano Zampini <[email protected]> ha scritto:

Are you trying to parallelize a 3 equations system? Or you just use your SIR code to experiment with TS?

...

On Mar 31, 2021, at 5:18 PM, Matthew Knepley <[email protected] <mailto:[email protected]>> wrote:

On Wed, Mar 31, 2021 at 10:15 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Thank you for your advices. I wrote what seems to me a very basic code, but I got this error when I run it with more than 1 processor: Clearly the result 299. is wrong but I do not understand what am doing wrong. With 1 processor it works fine.

My guess is that you do VecGetArray() and index the array using global indices rather than local indices, because there memory corruption with a Vec array.

Thanks,

Matt

steps 150, ftime 15. Vec Object: 2 MPI processes type: mpi Process [0] 16.5613 2.91405 Process [1] 299. [0]PETSC ERROR: PetscTrFreeDefault() called from VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [0]PETSC ERROR: Block [id=0(16)] at address 0x15812a0 is corrupted (probably write past end of array) [0]PETSC ERROR: Block allocated in VecCreate_MPI_Private() line 514 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pbvec.c [0]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [0]PETSC ERROR: Memory corruption: https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind <https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind> [0]PETSC ERROR: Corrupted memory [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Wed Mar 31 16:05:22 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 PetscTrFreeDefault() line 310 in /home/fbrarda/petsc/src/sys/memory/mtr.c [0]PETSC ERROR: #2 VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [0]PETSC ERROR: #3 VecDestroy() line 396 in /home/fbrarda/petsc/src/vec/vec/interface/vector.c [0]PETSC ERROR: #4 SNESLineSearchReset() line 284 in /home/fbrarda/petsc/src/snes/linesearch/interface/linesearch.c [0]PETSC ERROR: #5 SNESReset() line 3229 in /home/fbrarda/petsc/src/snes/interface/snes.c [0]PETSC ERROR: #6 TSReset() line 2800 in /home/fbrarda/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #7 TSDestroy() line 2856 in /home/fbrarda/petsc/src/ts/interface/ts.c [0]PETSC ERROR: #8 main() line 256 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries [0]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] <mailto:[email protected]>---------- application called MPI_Abort(MPI_COMM_SELF, 256001) - process 0 [1]PETSC ERROR: PetscTrFreeDefault() called from VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [1]PETSC ERROR: Block [id=0(16)] at address 0xbd9520 is corrupted (probably write past end of array) [1]PETSC ERROR: Block allocated in VecCreate_MPI_Private() line 514 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pbvec.c [1]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [1]PETSC ERROR: Memory corruption: https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind <https://www.mcs.anl.gov/petsc/documentation/installation.html#valgrind> [1]PETSC ERROR: Corrupted memory [1]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [1]PETSC ERROR: Petsc Release Version 3.14.4, unknown [1]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Wed Mar 31 16:05:22 2021 [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [1]PETSC ERROR: #1 PetscTrFreeDefault() line 310 in /home/fbrarda/petsc/src/sys/memory/mtr.c [1]PETSC ERROR: #2 VecDestroy_MPI() line 21 in /home/fbrarda/petsc/src/vec/vec/impls/mpi/pdvec.c [1]PETSC ERROR: #3 VecDestroy() line 396 in /home/fbrarda/petsc/src/vec/vec/interface/vector.c [1]PETSC ERROR: #4 TSReset() line 2806 in /home/fbrarda/petsc/src/ts/interface/ts.c [1]PETSC ERROR: #5 TSDestroy() line 2856 in /home/fbrarda/petsc/src/ts/interface/ts.c [1]PETSC ERROR: #6 main() line 256 in par_sir_model.c [1]PETSC ERROR: No PETSc Option Table entries [1]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] <mailto:[email protected]>---------- application called MPI_Abort(MPI_COMM_SELF, 256001) - process 0

...

Il giorno 31 mar 2021, alle ore 12:06, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Wed, Mar 31, 2021 at 3:54 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

You do not need one. We use it in examples because it makes it easy to create the data.

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Since you are producing all the Jacobian values, it is whatever is easier in your code I think.

THanks,

Matt

Best, Francesco

-- 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/ <https://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/>

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users <[email protected]> wrote:

...

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation. Barry

Hong (Mr.)

...

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda <[email protected]> wrote:

Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs. Thanks, Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/>

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda <[email protected]> wrote:

I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase". Matt

Thanks Francesco

Even though it will not scale and will deliver slower performance it is

...

completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda < [email protected]> wrote:

...

Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

...

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/>

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda < [email protected]> ha scritto:

In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda < [email protected]> wrote:

...

I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

...

Thanks Francesco

Even though it will not scale and will deliver slower performance it is

...

completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda < [email protected]> wrote:

...

Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

...

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/>

-- Stefano

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda <[email protected]> wrote:

Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan. Thanks, Matt

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini < [email protected]> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda < [email protected]> ha scritto:

...

In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda < [email protected]> wrote:

...

I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

...

Thanks Francesco

Even though it will not scale and will deliver slower performance it is

...

completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley < [email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda < [email protected]> wrote:

...

Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

...

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda < [email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/>

-- Stefano

-- 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/>

There are many examples that uses 1d DMDA in time dependent problems within the source tree. You can take a look at them first. Let us know if none of them do what you need. [szampini@localhost petsc]$ git grep DMDACreate1 src/ts/ src/ts/tests/ex12.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,appctx.m,1,1,NULL,&appctx.da);CHKERRQ(ierr); src/ts/tests/ex22.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,-n,1,1,0,&da);CHKERRQ(ierr); src/ts/tests/ex25.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,2,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/advection-diffusion-reaction/ex3.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/advection-diffusion-reaction/ex4.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE,8,2,1,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/advection-diffusion-reaction/ex6.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/ex10.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_NONE,M,1,0,lx,&da);CHKERRQ(ierr); src/ts/tutorials/ex10.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_NONE,20,sizeof(RDNode)/sizeof(PetscScalar),1,NULL,&rd->da);CHKERRQ(ierr); src/ts/tutorials/ex17.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,1,1,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/ex2.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,appctx.m,1,1,NULL,&appctx.da);CHKERRQ(ierr); src/ts/tutorials/ex21.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,appctx.m,1,1,NULL,&appctx.da);CHKERRQ(ierr); src/ts/tutorials/ex22.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,2,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/ex22f.F: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2, & src/ts/tutorials/ex22f_mf.F90: call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,im11,i2,i2,PETSC_NULL_INTEGER,da,ierr);CHKERRA(ierr) src/ts/tutorials/ex25.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,11,2,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/ex34.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 3, 1, NULL, &dm);CHKERRQ(ierr); src/ts/tutorials/ex4.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,m,1,1,NULL,&appctx.da);CHKERRQ(ierr); src/ts/tutorials/ex50.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,lenglob,1,1,NULL,&appctx.da);CHKERRQ(ierr); src/ts/tutorials/ex9.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/extchemfield.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,10,user.Nspec+1,1,NULL,& user.dm);CHKERRQ(ierr); src/ts/tutorials/multirate/ex5.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/multirate/ex6.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/multirate/ex7.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics.dof,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/multirate/ex8.c: ierr = DMDACreate1d(comm,DM_BOUNDARY_PERIODIC,50,ctx.physics2.dof,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/network/wash/pipeInterface.c: ierr = DMDACreate1d(PETSC_COMM_SELF, DM_BOUNDARY_GHOSTED, pipe->nnodes, 2, 1, NULL, &pipe->da);CHKERRQ(ierr); src/ts/tutorials/phasefield/biharmonic.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/phasefield/biharmonic2.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,2,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/phasefield/biharmonic3.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,2,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/phasefield/heat.c: ierr = DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, 10,1,2,NULL,&da);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9bus.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9bus.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busadj.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busadj.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busopt.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busopt.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busoptfd.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_gen,1,1,NULL,&user.dmgen);CHKERRQ(ierr); src/ts/tutorials/power_grid/stability_9bus/ex9busoptfd.c: ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,user.neqs_net,1,1,NULL,&user.dmnet);CHKERRQ(ierr); Il giorno gio 29 apr 2021 alle ore 10:57 Francesco Brarda < [email protected]> ha scritto:

Hi, The plan is actually to move to a SIR model also with the space. I understand that doing a SIR model in parallel will not bring any benefits, but I have been asked to do it as part of a project I am involved in.

I defined the DM as follows ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);

I am not sure whether I correctly understood this command properly. The vector should have 3 components (S, I, R) and 3 DOF as it is defined only when the three coordinates have been set. Then I create a global vector X. When I set the initial conditions as below

static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar *x; DM da;

PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr);

/* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);

x[0] = appctx->N - appctx->p[2]; x[1] = appctx->p[2]; x[2] = 0.0;

ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr); PetscFunctionReturn(0); }

I have the error:

[0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c [0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c [0]PETSC ERROR: #3 main() line 295 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries

I would be very happy to receive any advices to fix the code. Best, Francesco

Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda < [email protected]> wrote: Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.

Thanks,

Matt

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini < [email protected]> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda < [email protected]> ha scritto: In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda < [email protected]> wrote: I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

Thanks Francesco

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda < [email protected]> wrote: Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/

-- 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/

-- Stefano

-- 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/

-- Stefano

I am using DM to share the data when I define the functions outside the main (RHSFunction, InitialConditions and so on). In the sequential model, following the available examples, I used VecGetArray. As a replacement for a parallel model, I thought DM array would be a good alternative. Would you recommend using another tool? Best, Francesco

Il giorno 29 apr 2021, alle ore 15:27, Zhang, Hong <[email protected]> ha scritto:

Before you have a space-dependent SIR model, it does not make much sense to use DM and even consider parallelizing your code.

...

On Apr 29, 2021, at 2:57 AM, Francesco Brarda <[email protected] <mailto:[email protected]>> wrote:

Hi, The plan is actually to move to a SIR model also with the space. I understand that doing a SIR model in parallel will not bring any benefits, but I have been asked to do it as part of a project I am involved in.

I defined the DM as follows ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);

The size of DMDA is usually associated with the size of your computational domain, e.g. number of nodes in space. In your case, you have 3 degrees of freedom on each node, but you have only one node since the model has not been extended to the space yet.

Hong (Mr.)

...

I am not sure whether I correctly understood this command properly. The vector should have 3 components (S, I, R) and 3 DOF as it is defined only when the three coordinates have been set. Then I create a global vector X. When I set the initial conditions as below

static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar *x; DM da;

PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr);

/* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);

x[0] = appctx->N - appctx->p[2]; x[1] = appctx->p[2]; x[2] = 0.0;

ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr); PetscFunctionReturn(0); }

I have the error:

[0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c [0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c [0]PETSC ERROR: #3 main() line 295 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries

I would be very happy to receive any advices to fix the code. Best, Francesco

...

Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.

Thanks,

Matt

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini <[email protected] <mailto:[email protected]>> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda <[email protected] <mailto:[email protected]>> ha scritto: In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

Thanks Francesco

...

...

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

...

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected] <mailto:[email protected]>> ha scritto:

...

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users <[email protected] <mailto:[email protected]>> wrote:

...

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected] <mailto:[email protected]>> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Hong (Mr.)

...

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/ <https://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/ <https://www.cse.buffalo.edu/~knepley/>

-- Stefano

-- 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/ <https://www.cse.buffalo.edu/~knepley/>

[0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell Are you calling TSSetDM() to supply your created DMDA to the TS? Based on the error message you are not, it is using a default shell DM, which is what TS does if you do not provide them with a DM. You need to call TSSetDM() after you create the TS and the DMDA.

On Apr 29, 2021, at 2:57 AM, Francesco Brarda <[email protected]> wrote:

I defined the DM as follows ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);

If you truly want one "spatial" point then you would want to use 1,3,0. This says 1 location in space, three degrees of freedom at that point and 0 ghost values (since there is only one spatial point there can be no ghost spatial values). BUT DMDA ALWAYS puts all degree of freedom at a point on the same process so this will not give you parallelism. All 3 dof will be on the same MPI rank. For what you want to do you can use 3,1,2. This says three "spatial" points, 1 dof at each "spatial" point and 2 ghost points. In your case "spatial" does not mean spatial in space it is just three abstract points. The global vectors (from DMCreate or Get GlobalVector) will have 1 value on each rank. The local vector from DMCreate or Get LocalVector) will have three values on each rank. Your initial conditions code would be something like if (rank == 0) {

x[0] = appctx->N - appctx->p[2]; } else if (rank == 1) { x[1] = appctx->p[2]; } else { x[2] = 0.0;

} Your TSSetRHSFunction() would make a call to DMGetLocalVector(...&localX), do a DMGlobalToLocalBegin/End() from the input global X to localX, you would call DMDAVecGetArray(...,&xarray) on the localX and access all three values in xarray. The resulting computation of f the output vector would be something like if (rank == 0) {

farray[0] = your code that can use xarray[0], xarray[1], xarray[2] } else if (rank == 1) { farray[1] = your code that can use xarray[0], xarray[1], xarray[2] } else { farray[2] = your code that can use xarray[0], xarray[1], xarray[2]

} There are many examples of this pattern in the example tutorials. When you implement a code with a spatial distribution you would use a dof of 3 at each point and not parallelize over the dof at each point. Likely you want to use DMNETWORK to manage the spatial distribution since it has a simple API and allows any number of different number of neighbors for each point. DMDA would not make sense for true spatial distribution except in some truly trivial neighbor configurations. Barry

I am not sure whether I correctly understood this command properly. The vector should have 3 components (S, I, R) and 3 DOF as it is defined only when the three coordinates have been set. Then I create a global vector X. When I set the initial conditions as below

static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar *x; DM da;

PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr);

/* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);

x[0] = appctx->N - appctx->p[2]; x[1] = appctx->p[2]; x[2] = 0.0;

ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr); PetscFunctionReturn(0); }

I have the error:

[0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c [0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c [0]PETSC ERROR: #3 main() line 295 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries

I would be very happy to receive any advices to fix the code. Best, Francesco

...

Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.

Thanks,

Matt

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini <[email protected]> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda <[email protected]> ha scritto: In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda <[email protected]> wrote: I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

Thanks Francesco

...

...

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda <[email protected]> wrote: Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

...

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

...

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users <[email protected]> wrote:

...

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Hong (Mr.)

...

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/

-- 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/

-- Stefano

-- 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/

Using DMDA to get automatic parallelization of a 0-D code is nonsense. Ghost points are meant for spatially distributed data, not for 0D. If you really want to use DM, you should use DMRedundant and call DMGlobalToLocal to go to your distributed to local full representation (it calls MPI_Bcast internally). My advice is to read the many examples we have. Il giorno mar 4 mag 2021 alle ore 10:54 Francesco Brarda < [email protected]> ha scritto:

Thank you very much everyone. I do have one more question.

For what you want to do you can use 3,1,2. This says three "spatial" points, 1 dof at each "spatial" point and 2 ghost points. In your case "spatial" does not mean spatial in space it is just three abstract points.

In this case, since I have 2 ghost points, should I change the DMBoundaryType? For one process this works, but with 2 or 3 procs it doesn’t. The error I have is the following:

Solving a non-linear TS problem, number of processors = 2 [1]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [1]PETSC ERROR: Argument out of range [1]PETSC ERROR: Local x-width of domain x 1 is smaller than stencil width s 2 [1]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [1]PETSC ERROR: Petsc Release Version 3.14.4, unknown [1]PETSC ERROR: ./test_ic on a arch-debug named srvulx13 by fbrarda Mon May 3 17:46:37 2021 [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [1]PETSC ERROR: #1 DMSetUp_DA_1D() line 199 in /home/fbrarda/petsc/src/dm/impls/da/da1.c [1]PETSC ERROR: #2 DMSetUp_DA() line 20 in /home/fbrarda/petsc/src/dm/impls/da/dareg.c [1]PETSC ERROR: #3 DMSetUp() line 787 in /home/fbrarda/petsc/src/dm/interface/dm.c [1]PETSC ERROR: #4 main() line 232 in test_ic.c [1]PETSC ERROR: No PETSc Option Table entries [1]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected]

Il giorno 2 mag 2021, alle ore 18:54, Barry Smith <[email protected]> ha scritto:

[0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell

Are you calling TSSetDM() to supply your created DMDA to the TS? Based on the error message you are not, it is using a default shell DM, which is what TS does if you do not provide them with a DM. You need to call TSSetDM() after you create the TS and the DMDA.

On Apr 29, 2021, at 2:57 AM, Francesco Brarda <[email protected]> wrote:

I defined the DM as follows ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);

If you truly want one "spatial" point then you would want to use 1,3,0. This says 1 location in space, three degrees of freedom at that point and 0 ghost values (since there is only one spatial point there can be no ghost spatial values).

BUT DMDA ALWAYS puts all degree of freedom at a point on the same process so this will not give you parallelism. All 3 dof will be on the same MPI rank.

For what you want to do you can use 3,1,2. This says three "spatial" points, 1 dof at each "spatial" point and 2 ghost points. In your case "spatial" does not mean spatial in space it is just three abstract points.

The global vectors (from DMCreate or Get GlobalVector) will have 1 value on each rank. The local vector from DMCreate or Get LocalVector) will have three values on each rank. Your initial conditions code would be something like

if (rank == 0) {

x[0] = appctx->N - appctx->p[2];

} else if (rank == 1) {

x[1] = appctx->p[2];

} else {

x[2] = 0.0;

}

Your TSSetRHSFunction() would make a call to DMGetLocalVector(...&localX), do a DMGlobalToLocalBegin/End() from the input global X to localX, you would call DMDAVecGetArray(...,&xarray) on the localX and access all three values in xarray. The resulting computation of f the output vector would be something like

if (rank == 0) {

farray[0] = your code that can use xarray[0], xarray[1], xarray[2]

} else if (rank == 1) {

farray[1] = your code that can use xarray[0], xarray[1], xarray[2]

} else {

farray[2] = your code that can use xarray[0], xarray[1], xarray[2]

}

There are many examples of this pattern in the example tutorials.

When you implement a code with a spatial distribution you would use a dof of 3 at each point and not parallelize over the dof at each point. Likely you want to use DMNETWORK to manage the spatial distribution since it has a simple API and allows any number of different number of neighbors for each point. DMDA would not make sense for true spatial distribution except in some truly trivial neighbor configurations.

Barry

I am not sure whether I correctly understood this command properly. The vector should have 3 components (S, I, R) and 3 DOF as it is defined only when the three coordinates have been set. Then I create a global vector X. When I set the initial conditions as below

static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar *x; DM da;

PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr);

/* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);

x[0] = appctx->N - appctx->p[2]; x[1] = appctx->p[2]; x[2] = 0.0;

ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr); PetscFunctionReturn(0); }

I have the error:

[0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c [0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c [0]PETSC ERROR: #3 main() line 295 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries

I would be very happy to receive any advices to fix the code. Best, Francesco

Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda < [email protected]> wrote: Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.

Thanks,

Matt

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini < [email protected]> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda < [email protected]> ha scritto: In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda < [email protected]> wrote: I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

Thanks Francesco

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected]> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda < [email protected]> wrote: Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected]> ha scritto:

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users < [email protected]> wrote:

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected]> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

Hong (Mr.)

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/

-- 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/

-- Stefano

-- 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/

-- Stefano

Many thanks, this helped a lot. Now it finally works. Best, Francesco

Il giorno 4 mag 2021, alle ore 17:40, Barry Smith <[email protected]> ha scritto:

...

On May 4, 2021, at 2:53 AM, Francesco Brarda <[email protected] <mailto:[email protected]>> wrote:

Thank you very much everyone. I do have one more question.

...

For what you want to do you can use 3,1,2. This says three "spatial" points, 1 dof at each "spatial" point and 2 ghost points. In your case "spatial" does not mean spatial in space it is just three abstract points. In this case, since I have 2 ghost points, should I change the DMBoundaryType? For one process this works, but with 2 or 3 procs it doesn’t. The error I have is the following:

Solving a non-linear TS problem, number of processors = 2 [1]PETSC ERROR: --------------------- Error Message -------------------------------------------------------------- [1]PETSC ERROR: Argument out of range [1]PETSC ERROR: Local x-width of domain x 1 is smaller than stencil width s 2

I had forgotten this. It is a limitation of the DMDA implementation that one cannot require data from two ranks away from the current rank.

But yes perhaps you can "cheat" by using DM_BOUNDARY_PERIODIC and a stencil width of 1. For each rank it will bring the value from the previous and next rank and due to the periodicity it will thus bring the 3rd value to the first rank and the first value to the third rank.

Barry

...

[1]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [1]PETSC ERROR: Petsc Release Version 3.14.4, unknown [1]PETSC ERROR: ./test_ic on a arch-debug named srvulx13 by fbrarda Mon May 3 17:46:37 2021 [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [1]PETSC ERROR: #1 DMSetUp_DA_1D() line 199 in /home/fbrarda/petsc/src/dm/impls/da/da1.c [1]PETSC ERROR: #2 DMSetUp_DA() line 20 in /home/fbrarda/petsc/src/dm/impls/da/dareg.c [1]PETSC ERROR: #3 DMSetUp() line 787 in /home/fbrarda/petsc/src/dm/interface/dm.c [1]PETSC ERROR: #4 main() line 232 in test_ic.c [1]PETSC ERROR: No PETSc Option Table entries [1]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected] <mailto:[email protected]>----------

...

Il giorno 2 mag 2021, alle ore 18:54, Barry Smith <[email protected] <mailto:[email protected]>> ha scritto:

[0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell

Are you calling TSSetDM() to supply your created DMDA to the TS? Based on the error message you are not, it is using a default shell DM, which is what TS does if you do not provide them with a DM. You need to call TSSetDM() after you create the TS and the DMDA.

...

On Apr 29, 2021, at 2:57 AM, Francesco Brarda <[email protected] <mailto:[email protected]>> wrote:

I defined the DM as follows ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);

If you truly want one "spatial" point then you would want to use 1,3,0. This says 1 location in space, three degrees of freedom at that point and 0 ghost values (since there is only one spatial point there can be no ghost spatial values).

BUT DMDA ALWAYS puts all degree of freedom at a point on the same process so this will not give you parallelism. All 3 dof will be on the same MPI rank.

For what you want to do you can use 3,1,2. This says three "spatial" points, 1 dof at each "spatial" point and 2 ghost points. In your case "spatial" does not mean spatial in space it is just three abstract points.

The global vectors (from DMCreate or Get GlobalVector) will have 1 value on each rank. The local vector from DMCreate or Get LocalVector) will have three values on each rank. Your initial conditions code would be something like

if (rank == 0) {

...

x[0] = appctx->N - appctx->p[2]; } else if (rank == 1) { x[1] = appctx->p[2]; } else { x[2] = 0.0; }

Your TSSetRHSFunction() would make a call to DMGetLocalVector(...&localX), do a DMGlobalToLocalBegin/End() from the input global X to localX, you would call DMDAVecGetArray(...,&xarray) on the localX and access all three values in xarray. The resulting computation of f the output vector would be something like

if (rank == 0) {

...

farray[0] = your code that can use xarray[0], xarray[1], xarray[2] } else if (rank == 1) { farray[1] = your code that can use xarray[0], xarray[1], xarray[2] } else { farray[2] = your code that can use xarray[0], xarray[1], xarray[2] }

There are many examples of this pattern in the example tutorials.

When you implement a code with a spatial distribution you would use a dof of 3 at each point and not parallelize over the dof at each point. Likely you want to use DMNETWORK to manage the spatial distribution since it has a simple API and allows any number of different number of neighbors for each point. DMDA would not make sense for true spatial distribution except in some truly trivial neighbor configurations.

Barry

...

I am not sure whether I correctly understood this command properly. The vector should have 3 components (S, I, R) and 3 DOF as it is defined only when the three coordinates have been set. Then I create a global vector X. When I set the initial conditions as below

static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar *x; DM da;

PetscFunctionBeginUser; ierr = TSGetDM(ts,&da);CHKERRQ(ierr);

/* Get pointers to vector data */ ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);

x[0] = appctx->N - appctx->p[2]; x[1] = appctx->p[2]; x[2] = 0.0;

ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr); PetscFunctionReturn(0); }

I have the error:

[0]PETSC ERROR: Invalid argument [0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell [0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html <https://www.mcs.anl.gov/petsc/documentation/faq.html> for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.14.4, unknown [0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021 [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --with-openblas-dir=/opt/packages/openblas/0.2.13-gcc --download-mpich PETSC_ARCH=arch-debug [0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c [0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c [0]PETSC ERROR: #3 main() line 295 in par_sir_model.c [0]PETSC ERROR: No PETSc Option Table entries

I would be very happy to receive any advices to fix the code. Best, Francesco

...

Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized. Do you think it is better to put the 3 equations into a 1d Distributed Array with 3 dofs and run the job with multiple procs regardless of how many equations I have? Is it possible?

If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.

Thanks,

Matt

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini <[email protected] <mailto:[email protected]>> ha scritto:

It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly. Is this work you are doing used to understand PETSc parallelization strategy? if so, there are multiple examples in the sourcetree that you can look at to populate matrices and vectors in parallel

Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda <[email protected] <mailto:[email protected]>> ha scritto: In principle the entire code was for 1 proc only. The functions were built with VecGetArray(). While adapting the code for multiple procs I thought using VecGetOwnershipRange was a possible way to allocate the equations in the vector using multiple procs. What do you think, please?

Thank you, Francesco

...

Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: I was trying to follow Barry's advice some time ago, but I guess that's not the way he meant it. How should I refer to the values contained in x? With Distributed Arrays?

That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".

Matt

Thanks Francesco

...

...

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:

On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda <[email protected] <mailto:[email protected]>> wrote: Hi! I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment. I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:

I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only 1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.

Thanks,

Matt

static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscScalar f;//, *x_localptr; const PetscScalar *x; PetscInt mybase;

PetscFunctionBeginUser; ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 1) { f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } if (mybase == 2) { f = (PetscScalar) (appctx->p2*x[1]); ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr); ierr = VecAssemblyBegin(F);CHKERRQ(ierr); ierr = VecAssemblyEnd(F);CHKERRQ(ierr); PetscFunctionReturn(0); }

Whilst for the Jacobian I did:

static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx) { PetscErrorCode ierr; AppCtx *appctx = (AppCtx*) ctx; PetscInt mybase, rowcol[] = {0,1,2}; const PetscScalar *x;

PetscFunctionBeginUser; ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr); ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr); if (mybase == 0) { const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 1) { const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } if (mybase == 2) { const PetscScalar J[] = {0, - appctx->p2, a}; ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr); } ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); if (A != B) { ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); } PetscFunctionReturn(0); }

This code does not provide the correct result, that is, the solution is the initial condition, either using implicit or explicit methods. Is the way I defined these objects wrong? How can I fix it? I also tried to print the Jacobian with the following commands but it does not work (blank rows and error message). How should I print the Jacobian?

ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr); ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);

I would very much appreciate any kind of help or advice. Best, Francesco

...

Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <[email protected] <mailto:[email protected]>> ha scritto:

...

On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users <[email protected] <mailto:[email protected]>> wrote:

...

On Mar 31, 2021, at 2:53 AM, Francesco Brarda <[email protected] <mailto:[email protected]>> wrote:

Hi everyone!

I am trying to solve a system of 3 ODEs (a basic SIR model) with TS. Sequentially works pretty well, but I need to switch it into a parallel version. I started working with TS not very long time ago, there are few questions I’d like to share with you and if you have any advices I’d be happy to hear. First of all, do I need to use a DM object even if the model is only time dependent? All the examples I found were using that object for the other variable when solving PDEs.

Are you considering SIR on a spatial domain? If so, you can parallelize your model in the spatial domain using DM. Splitting the three variables in the ODE among processors would not scale.

Even though it will not scale and will deliver slower performance it is completely possible for you to solve the 3 variable problem using 3 MPI ranks. Or 10 mpi ranks. You would just create vectors/matrices with 1 degree of freedom for the first three ranks and no degrees of freedom for the later ranks. During your function evaluation (and Jacobian evaluation) for TS you will need to set up the appropriate communication to get the values you need on each rank to evaluate the parts of the function evaluation needed by that rank. This is true for parallelizing any computation.

Barry

...

Hong (Mr.)

...

When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?

Best, Francesco

-- 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/ <https://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/ <https://www.cse.buffalo.edu/~knepley/>

-- Stefano

-- 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/ <https://www.cse.buffalo.edu/~knepley/>