<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">I am using DM to share the data when I define the functions outside the main (RHSFunction, InitialConditions and so on).<div class="">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?</div><div class=""><br class=""></div><div class="">Best,</div><div class="">Francesco</div><div class=""><br class=""><div><br class=""><blockquote type="cite" class=""><div class="">Il giorno 29 apr 2021, alle ore 15:27, Zhang, Hong <<a href="mailto:hongzhang@anl.gov" class="">hongzhang@anl.gov</a>> ha scritto:</div><br class="Apple-interchange-newline"><div class="">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" class="">
<div style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">
Before you have a space-dependent SIR model, it does not make much sense to use DM and even consider parallelizing your code. <br class="">
<div class=""><br class="">
<blockquote type="cite" class="">
<div class="">On Apr 29, 2021, at 2:57 AM, Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> wrote:</div>
<br class="Apple-interchange-newline">
<div class="">
<div style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">
Hi,<br class="">
The plan is actually to move to a SIR model also with the space.
<div class="">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.</div>
<div class=""><br class="">
I defined the DM as follows<br class="">
<font face="Menlo" class="">ierr = DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,3,3,3,NULL,&da);CHKERRQ(ierr);</font><br class="">
</div>
</div>
</div>
</blockquote>
<div class=""><br class="">
</div>
<div class="">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.</div>
<div class=""><br class="">
</div>
<div class="">Hong (Mr.)</div>
<br class="">
<blockquote type="cite" class="">
<div class="">
<div style="word-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;" class="">
<div class=""><br class="">
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.<br class="">
Then I create a global vector X. When I set the initial conditions as below
<div class=""><br class="">
<font face="Menlo" class="">static PetscErrorCode InitialConditions(TS ts,Vec X, void *ctx)<br class="">
{ <br class="">
PetscErrorCode ierr;<br class="">
AppCtx *appctx = (AppCtx*) ctx;<br class="">
PetscScalar *x;<br class="">
DM da;<br class="">
<br class="">
PetscFunctionBeginUser;<br class="">
ierr = TSGetDM(ts,&da);CHKERRQ(ierr);<br class="">
<br class="">
/* Get pointers to vector data */<br class="">
ierr = DMDAVecGetArray(da,X,(void*)&x);CHKERRQ(ierr);<br class="">
<br class="">
x[0] = appctx->N - appctx->p[2];<br class="">
x[1] = appctx->p[2];<br class="">
x[2] = 0.0;<br class="">
<br class="">
ierr = DMDAVecRestoreArray(da,X,(void*)&x);CHKERRQ(ierr);<br class="">
PetscFunctionReturn(0);<br class="">
}<br class="">
</font><br class="">
I have the error:<br class="">
<br class="">
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures;" class="">[0]PETSC ERROR: Invalid argument</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: Wrong subtype object:Parameter # 1 must have implementation da it is shell</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: See
<a href="https://www.mcs.anl.gov/petsc/documentation/faq.html" class="">https://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: Petsc Release Version 3.14.4, unknown </span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: ./par_sir_model on a arch-debug named srvulx13 by fbrarda Thu Apr 29 09:36:17 2021</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[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</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: #1 DMDAVecGetArray() line 48 in /home/fbrarda/petsc/src/dm/impls/da/dagetarray.c</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: #2 InitialConditions() line 175 in par_sir_model.c</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: #3 main() line 295 in par_sir_model.c</span></div>
<div style="margin: 0px; font-stretch: normal; line-height: normal; font-family: Menlo;" class="">
<span style="font-variant-ligatures: no-common-ligatures" class="">[0]PETSC ERROR: No PETSc Option Table entries</span></div>
<div class=""><br class="">
</div>
<div class="">I would be very happy to receive any advices to fix the code.</div>
<div class="">Best,</div>
<div class="">Francesco</div>
<br class="">
<blockquote type="cite" class="">Il giorno 20 apr 2021, alle ore 21:35, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> ha scritto:<br class="">
<br class="">
On Tue, Apr 20, 2021 at 1:17 PM Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> wrote:<br class="">
Thank you for the advices, I would just like to convert the code I already have to see what might happen once parallelized.<br class="">
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? <br class="">
<br class="">
If you plan in the end to use a structured grid, this is a great plan. If not, this is not a good plan.<br class="">
<br class="">
Thanks,<br class="">
<br class="">
Matt<br class="">
<br class="">
Thank you,<br class="">
Francesco<br class="">
<br class="">
<blockquote type="cite" class="">Il giorno 20 apr 2021, alle ore 17:57, Stefano Zampini <<a href="mailto:stefano.zampini@gmail.com" class="">stefano.zampini@gmail.com</a>> ha scritto:<br class="">
<br class="">
It does not make sense to parallelize to 1 equation per process, unless that single equation per process is super super super costly.<br class="">
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<br class="">
<br class="">
Il giorno mar 20 apr 2021 alle ore 17:52 Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> ha scritto:<br class="">
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? <br class="">
<br class="">
Thank you,<br class="">
Francesco<br class="">
<br class="">
<blockquote type="cite" class="">Il giorno 20 apr 2021, alle ore 16:43, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> ha scritto:<br class="">
<br class="">
On Tue, Apr 20, 2021 at 10:41 AM Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> wrote:<br class="">
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?<br class="">
<br class="">
That is how you get values from x. However, I cannot understand at all what you are doing with "mybase".<br class="">
<br class="">
Matt<br class="">
<br class="">
Thanks <br class="">
Francesco<br class="">
<br class="">
<blockquote type="cite" class="">
<blockquote type="cite" class=""> 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.<br class="">
<br class="">
Barry<br class="">
</blockquote>
</blockquote>
<br class="">
<br class="">
<br class="">
<blockquote type="cite" class="">Il giorno 20 apr 2021, alle ore 15:40, Matthew Knepley <<a href="mailto:knepley@gmail.com" class="">knepley@gmail.com</a>> ha scritto:<br class="">
<br class="">
On Tue, Apr 20, 2021 at 9:36 AM Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> wrote:<br class="">
Hi!<br class="">
I tried to implement the SIR model taking into account the fact that I will only use 3 MPI ranks at this moment.<br class="">
I built vectors and matrices following the examples already available. In particular, I defined the functions required similarly (RHSFunction, IFunction, IJacobian), as follows:<br class="">
<br class="">
I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only<br class="">
1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.<br class="">
<br class="">
Thanks,<br class="">
<br class="">
Matt<br class="">
<br class="">
static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ctx)<br class="">
{ <br class="">
PetscErrorCode ierr;<br class="">
AppCtx *appctx = (AppCtx*) ctx;<br class="">
PetscScalar f;//, *x_localptr; <br class="">
const PetscScalar *x;<br class="">
PetscInt mybase;<br class="">
<br class="">
PetscFunctionBeginUser;<br class="">
ierr = VecGetOwnershipRange(X,&mybase,NULL);CHKERRQ(ierr);<br class="">
ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);<br class="">
if (mybase == 0) {<br class="">
f = (PetscScalar) (-appctx->p1*x[0]*x[1]/appctx->N);<br class="">
ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES);<br class="">
}<br class="">
if (mybase == 1) {<br class="">
f = (PetscScalar) (appctx->p1*x[0]*x[1]/appctx->N-appctx->p2*x[1]);<br class="">
ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES);<br class="">
}<br class="">
if (mybase == 2) {<br class="">
f = (PetscScalar) (appctx->p2*x[1]);<br class="">
ierr = VecSetValues(F,1,&mybase,&f,INSERT_VALUES);<br class="">
}<br class="">
ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);<br class="">
ierr = VecAssemblyBegin(F);CHKERRQ(ierr);<br class="">
ierr = VecAssemblyEnd(F);CHKERRQ(ierr);<br class="">
PetscFunctionReturn(0);<br class="">
}<br class="">
<br class="">
<br class="">
Whilst for the Jacobian I did:<br class="">
<br class="">
<br class="">
static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)<br class="">
{ <br class="">
PetscErrorCode ierr;<br class="">
AppCtx *appctx = (AppCtx*) ctx;<br class="">
PetscInt mybase, rowcol[] = {0,1,2};<br class="">
const PetscScalar *x;<br class="">
<br class="">
PetscFunctionBeginUser;<br class="">
ierr = MatGetOwnershipRange(B,&mybase,NULL);CHKERRQ(ierr);<br class="">
ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);<br class="">
if (mybase == 0) {<br class="">
const PetscScalar J[] = {a + appctx->p1*x[1]/appctx->N, appctx->p1*x[0]/appctx->N, 0};<br class="">
ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr);<br class="">
}<br class="">
if (mybase == 1) {<br class="">
const PetscScalar J[] = {- appctx->p1*x[1]/appctx->N, a - appctx->p1*x[0]/appctx->N + appctx->p2, 0};<br class="">
ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr);<br class="">
}<br class="">
if (mybase == 2) {<br class="">
const PetscScalar J[] = {0, - appctx->p2, a};<br class="">
ierr = MatSetValues(B,1,&mybase,3,rowcol,J,INSERT_VALUES);CHKERRQ(ierr);<br class="">
}<br class="">
ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);<br class="">
<br class="">
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
if (A != B) {<br class="">
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
}<br class="">
PetscFunctionReturn(0);<br class="">
}<br class="">
<br class="">
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? <br class="">
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?<br class="">
<br class="">
ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL); CHKERRQ(ierr);<br class="">
ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class="">
ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); <br class="">
<br class="">
I would very much appreciate any kind of help or advice.<br class="">
Best,<br class="">
Francesco<br class="">
<br class="">
<blockquote type="cite" class="">Il giorno 2 apr 2021, alle ore 04:45, Barry Smith <<a href="mailto:bsmith@petsc.dev" class="">bsmith@petsc.dev</a>> ha scritto:<br class="">
<br class="">
<br class="">
<br class="">
<blockquote type="cite" style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none;" class="">
On Apr 1, 2021, at 9:17 PM, Zhang, Hong via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" class="">petsc-users@mcs.anl.gov</a>> wrote:<br class="">
<br class="">
<br class="">
<br class="">
<blockquote type="cite" class="">On Mar 31, 2021, at 2:53 AM, Francesco Brarda <<a href="mailto:brardafrancesco@gmail.com" class="">brardafrancesco@gmail.com</a>> wrote:<br class="">
<br class="">
Hi everyone!<br class="">
<br class="">
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. <br class="">
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.<br class="">
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.<br class="">
</blockquote>
<br class="">
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.<br class="">
</blockquote>
<br class="">
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.<br class="">
<br class="">
Barry<br class="">
<br class="">
<br class="">
<br class="">
<br class="">
<blockquote type="cite" style="font-family: Helvetica; font-size: 12px; font-style: normal; font-variant-caps: normal; font-weight: normal; letter-spacing: normal; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; word-spacing: 0px; text-decoration: none;" class="">
<br class="">
Hong (Mr.) <br class="">
<br class="">
<blockquote type="cite" class="">When I preallocate the space for the Jacobian matrix, is it better to decide the local or global space?<br class="">
<br class="">
Best,<br class="">
Francesco<br class="">
</blockquote>
</blockquote>
</blockquote>
<br class="">
<br class="">
<br class="">
-- <br class="">
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">
-- Norbert Wiener<br class="">
<br class="">
<a href="https://www.cse.buffalo.edu/~knepley/" class="">https://www.cse.buffalo.edu/~knepley/</a><br class="">
</blockquote>
<br class="">
<br class="">
<br class="">
-- <br class="">
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">
-- Norbert Wiener<br class="">
<br class="">
<a href="https://www.cse.buffalo.edu/~knepley/" class="">https://www.cse.buffalo.edu/~knepley/</a><br class="">
</blockquote>
<br class="">
<br class="">
<br class="">
-- <br class="">
Stefano<br class="">
</blockquote>
<br class="">
<br class="">
<br class="">
-- <br class="">
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br class="">
-- Norbert Wiener<br class="">
<br class="">
<a href="https://www.cse.buffalo.edu/~knepley/" class="">https://www.cse.buffalo.edu/~knepley/</a><br class="">
</blockquote>
<br class="">
</div>
</div>
</div>
</div>
</blockquote>
</div>
<br class="">
</div>
</div></blockquote></div><br class=""></div></body></html>