[petsc-users] a question about DAE Function

[Jin, Shuangshuang]

I see. So for my IFunction which is called by TSSetIFunction, I can leave it as what it is right now (I mean scatter x and define equations), but for the IJacobian function which is called by TSSetIJacobian, I should define the Jacobian matrix as MATELEMENTAL type and set up its values accordingly? 

My Jacobian matrix also use the values of the distributed array X. Generally, the benefit of MATELEMENTAL is to avoid the scattering of X to each processor inside the IJacobian function in each iteration by "using MatGetOwnershipIS() to determine which matrix entries are responsible for computing" and such?

Thanks,
Shuangshuang

-----Original Message-----
From: Jed Brown [mailto:five9a2 at gmail.com] On Behalf Of Jed Brown
Sent: Wednesday, June 05, 2013 12:02 PM
To: Jin, Shuangshuang
Cc: petsc-users at mcs.anl.gov
Subject: RE: [petsc-users] a question about DAE Function

"Jin, Shuangshuang" <Shuangshuang.Jin at pnnl.gov> writes:

> Hi, thanks for your message, but I still have a few questions here.
>
> 1. How to "reuse the VecScatter"?

Store it in your user context and reuse it the next time you need to scatter.

> 2. The data movement of O(N) is for scattering X vector to each 
> processor? And the computation of O(N^3) is for solving the DAE 
> equation in parallel? So "the parallelism will be useful if N is large 
> enough" as you indicated. Then does that mean my current 
> implementation of scattering X and defining f[] on different 
> processors to solve the DAE in parallel is the right way to do it?
> Then what's the benefit of using MatElemental?

That's fine to define the equations, but you still have to solve them somehow.  Most common methods do that using a matrix.

> 3. When I looked into the example of MatElemental, it seems to me that 
> they're all matrix operations. I cannot understand its connection with 
> my code, especially the IFUNCTION part. I only have several functions 
> definitions f[] using the data of a distributed array X. Don't know 
> where and how to use MatElemental here.

Are you assembling a Jacobian (TSSetIJacobian)?  If not, then Elemental will not help, but keep in mind that this will converge very slowly unless your system is well-conditioned.