[petsc-users] low dimensional sub-problem

[Jed Brown]

Gideon Simpson <gideon.simpson at gmail.com> writes:

> I’ve been continuing working on implementing a projection method problem, which, loosely, looks like the following.  Vec y contains the state vector for my system, y’ = f(y) which is solved with a TS, using, for now, rk4.
>
> I have added to this a TSPostStep which does root finding on the nonlinear problem
>
> f(lambda) = f(lambda; y, w) =  g(y + lambda * w) - g0

You might want to treat this as a DAE.

> y and w are vectors that are the size of the system (numbers of mesh points), and lambda is a scalar perturbation.  
>
> Right now, the snes function for solving this amounts to:
>
> PetscErrorCode projector(SNES snes, Vec lam, Vec f, void *ctx){
>
>   AppCtx * user = (AppCtx *) ctx;
>   const PetscScalar *lam_value;
>   PetscScalar *f_value;
>
>   VecGetArrayRead(lam, &lam_value);
>   VecGetArray(f, &f_value);
>   VecCopy(user->y, user->w);
>   VecAXPY(user->w, lam_value[0], user->gradMy);
>
>   f_value[0] = g(user->w) -user->g0;
>   VecRestoreArrayRead(lam, &lam_value);
>   VecRestoreArray(f, &f_value);
>   
>   return 0;
> }
>
> To get this all to work, I constructed the SNES and Vec lam with PETSC_COMM_SELF, effectively putting a copy of the nonlinear problem on each process.  Basically, the nonlinear problem is low dimensional, but it parametrically depends on the high dimensional, distributed, vectors y and w.
>
> The one thing that bothers me about this is:  
>
> 1.  It would seem that each process is is doing all of these vector operations, which is entirely redundant, as the only thing that really needs to be shared amongst the processes is the value of lambda.  Is that correct? 

Only the SNES part is redundant, but it's super cheap because it's a
scalar problem.  The Vec operations are using global (parallel,
distributed) vectors, so they are not redundant.  Of course it is
critical for the SNES on every process to be configured identically and
deterministic so that the processes don't diverge.  And g(user->w) must
return the same value on each process (it probably finishes with an
MPI_Allreduce or similar).

> 2.  Is there an obvious way around this redundancy?
>
> -gideon
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 832 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170225/c7e88aad/attachment.pgp>