[petsc-users] a question about DAE Function

[Jin, Shuangshuang]

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)
      }
}

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, but the problem is fully coupled.

  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?

  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?

Thanks,
Shuangshuang

   Matt

Where x is an distributed array, a and c are constant arrays, b and d are constants. n is the length of the global length of x.

The problem is, since x is distributed across all the processes, each process only owns a piece of items of x, when I do (1) as written above, no matter which processor f[i+2] is on, it has to grab all the x array data from other processors, which means I have to do a scattering of x array in advance to make the x[k] value available. This is quite an inefficient way given the scattering of vector would happen in each iteration when the IFuntion is called.

I am wondering if this is the only way I have to deal with this problem, or am I trapped in a wrong direction?

Your comments are highly appreciated!

Thanks,
Shuangshuang





--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130604/ef5c713c/attachment.html>