<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 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?<div class=""><br class=""><div class="">Thanks </div><div class="">Francesco</div><div class=""><br class=""></div><div class=""><blockquote type="cite" class=""><div dir="ltr" class=""><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div class="" style="overflow-wrap: break-word;"><div class=""><span 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></span></div></div></blockquote></div></div></blockquote><div class=""><br class=""></div><br class=""><div><br class=""><blockquote type="cite" class=""><div 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:</div><br class="Apple-interchange-newline"><div class=""><meta charset="UTF-8" class=""><div dir="ltr" style="caret-color: rgb(0, 0, 0); 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; -webkit-text-stroke-width: 0px; text-decoration: none;" class=""><div dir="ltr" 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=""></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="overflow-wrap: break-word;" 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=""></div></blockquote><div class=""><br class=""></div><div class="">I don't think this makes sense. You use "mybase" to distinguish between 3 procs, which would indicate that each procs has only</div><div class="">1 degree of freedom. However, you use x[1] on each proc, indicating it has at least 2 dofs.</div><div class=""><br class=""></div><div class=""> Thanks,</div><div class=""><br class=""></div><div class=""> Matt</div><div class=""> </div><blockquote class="gmail_quote" style="margin: 0px 0px 0px 0.8ex; border-left-width: 1px; border-left-style: solid; border-left-color: rgb(204, 204, 204); padding-left: 1ex;"><div style="overflow-wrap: break-word;" class=""><font face="Menlo" 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="">}</font><br class=""><br class=""><br class="">Whilst for the Jacobian I did:<div class=""><br class=""></div><div class=""><br class=""><font face="Menlo" 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="">}</font><div class=""><font face="Menlo" class=""><br class=""></font></div><div 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? </div><div 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?</div><div class=""><br class=""></div><font face="Menlo" class=""><span class="">ierr = TSGetIJacobian(ts,NULL,&K, NULL, NULL);</span></font> <span style="font-family: Menlo;" class="">CHKERRQ(ierr);</span><br style="font-family: Menlo;" class=""><font face="Menlo" class=""><span class="">ierr = MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class=""></span><span class="">ierr = MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);<br class=""></span></font><span class=""><font face="Menlo" class="">ierr = MatView(K,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</font> <br class=""><br class="">I would very much appreciate any kind of help or advice.</span></div><div class=""><span class="">Best,</span></div><div class=""><span 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" target="_blank" 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" target="_blank" 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" target="_blank" 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=""></span></div></div></blockquote></div><br clear="all" class=""><div class=""><br class=""></div>--<span class="Apple-converted-space"> </span><br class=""><div dir="ltr" class="gmail_signature"><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div class=""><div dir="ltr" class=""><div 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</div><div class=""><br class=""></div><div class=""><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank" class="">https://www.cse.buffalo.edu/~knepley/</a></div></div></div></div></div></div></div></div></div></blockquote></div><br class=""></div></div></body></html>