[petsc-users] a question about DAE Function

[Jed Brown]

Your message quoting is very peculiar.  Please use normal quoting in the
future.

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

> From: Matthew Knepley [mailto:knepley at gmail.com]
> Sent: Tuesday, June 04, 2013 12:32 PM
> To: Jin, Shuangshuang
> Cc: petsc-users at mcs.anl.gov
> Subject: Re: [petsc-users] a question about DAE Function
>
> On Tue, Jun 4, 2013 at 3:13 PM, Jin, Shuangshuang <Shuangshuang.Jin at pnnl.gov<mailto:Shuangshuang.Jin at pnnl.gov>> wrote:
> Hello, I have a question regarding how to set up the IFunction for TSSetIFunction. Please see the detailed illustration below:
>
> Inside my IFunction, I have one of the f equations defined as:
>
> if (i >= xstart && i < xstart+xlen) {
>       f[i+2] = sin(x[i])*x[i+2] + cos(x[i])*x[i+3];
>
>
> 1)      You should not be setting f[i+2]. You always set the residual local to your process.
> [Jin, Shuangshuang]  Sorry, I typed it wrong, the condition is "if (i+2 >= xstart && i+2 < xstart+xlen) {" instead.
>
>       for (k = 0; k < n; k++) {
>           f[i+2] += -(cos(x[k])*a[k])*b) - (sin(x[k])*c[k])*d);               (1)
>       }
> }

Showing something of the outer loop would have made the usage more
clear.

>
> 2) This in effect says that your Jacobian is completely dense, or that your problem is fully coupled. I see two
>     choices:
> [Jin, Shuangshuang] I have four sets of f functions: f[i], f[i+1],
>     f[i+2], and f[i+3]. (i=1..n), the Jacobian matrix for the whole
>     DAE problem is sparse, 

Your code above couples f[i+2] to x[k] for k=1..n, which is dense.

> but the problem is fully coupled.

"fully coupled" means "dense", unless you intended something different.
What exactly is sparse?

>   a) Use scalable dense linear algebra, like Elemental, and formulate your operations as linear algebra
> [Jin, Shuangshuang] Do you have a sample source code on the web for this?

There are a couple tutorials using MatElemental, but we don't seem to
agree about whether the problem formulation is correct, so don't do that
yet.

>   b) If this is a boundary integral equation, consider using the FMM technology
> [Jin, Shuangshuang] It's not.
>
> Either way, this is not a PDE, and so using that PETSc parallel paradigm is wrong.
> [Jin, Shuangshuang] This is a DAE, I set up the section 6.1.2 Solving Differential Algebraic Equations of the Manual, form IFUNCITON, IJacobain, and set up initial values for x. If I scatter x every time inside IFUNCTION, the problem is solved in parallel when I use multiple threads. I just felt that my way of scattering x is inefficient. I don't understand why using the PETSc parallel paradigm to solve the problem is wrong?

If your function is dense, then you need to gather all of x to compute a
residual.  That is what your code above implies, but you also said it
was sparse, so which is it?