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