Working between two boundaries

[tipara]

Hi.
I'm trying to use OpenCurrent for a school project, but I don't quite
know how to. Basically I need to solve Laplace's equation between the
boundaries of the grid, and a body in the middle of the grid, with
Neumann's conditions on both boundaries. For the first I assume there
is no difficulty using BC_NEUMANN, but I don't see how I can do the
same on the edges of the body. Is there any way to do it ?

Also, I'm not sure I've understood correctly how to use solvers : are
the correct steps :
1/ initialize the solver and grids
2/ fill the grid with the values of phi and (d phi/d t) I want to
enforce, set n{x,y,z} and h{x,y,z}
3/ set parameters and then solve, then read from the solver using
copy_all_data ?

Sorry if it was not the right place to ask, and thank you.

[Jonathan Cohen]

Yes, that is basically the correct way to use the solvers.

http://code.google.com/p/opencurrent/source/browse/src/tests/diffusion3dtest.cpp

is probably a good template to follow.  The run_tests() routine basically does what you say.

As for how to add internal Neumann conditions, there are a couple of ways.  For Neumann conditions, probably the best way is to write a modified version of the Sol_LaplacianCentered3DDevice that actually enforces the Neumann conditions while computing the Laplacian operator.  To figure out where the body is, you need to add a new grid:

  Grid3DDevice<int> _flags;

_flags would be 0 for grid cells that are not inside the body, 1 for grid cells that are inside.  (BTW, I'm assuming the Neumann condition is dphi/dn = 0)

Then inside the solve() routine, you call a kernel that looks like basically the same, but only adds each directional term to the Laplacian if the Neumann condition flag for that cell is not set:

template<typename T>
__global__ void Sol_LaplacianCentered3DDevice_stencil_with_boundary(
  T *phi, T *dphidt, const int *flags, ... (same as before)
{
  int i,j,k = ...; // same as before

  if (i < nx && j < ny && k < nz) {

    // calc phi indexing
    int idx = __mul24(i, xstride) + __mul24(j,ystride) + k;
   
    T phi_ijk = phi[idx];
    T dphi_ijk = dphidt[idx];

    T laplacian =
      flags[idx + 1] ? 0 : invhz2 * (phi[idx + 1] - phi_ijk) +
      flags[idx - 1] ? 0 : invhz2 * (phi[idx - 1] - phi_ijk) +
      flags[idx + ystride] ? 0 : invhy2 * (phi[idx + ystride] - phi_ijk) +
      flags[idx - ystride] ? 0 : invhy2 * (phi[idx - ystride ] - phi_ijk) +
      flags[idx + xstride] ? 0 : invhx2 * (phi[idx + xstride] - phi_ijk) +
      flags[idx - xstride] ? 0 : invhx2 * (phi[idx - xstride ] - phi_ijk);

    dphidt[idx] = dphi_ijk + coefficient * laplacian;
  }
}


Does this make sense?

[tipara]

Thank you, I've managed to get the modified code to run, but I was
puzzled by the results of the computation, until I (finally) realized
that the equation being solved was d phi/dt = coefficient*(laplacian
phi), not laplacian phi = 0... Is there any way to solve the latter
(the goal being to simulate a potential airflow around an object) ?

Anyway, thank you very much for your help.

[dpephd]

A linear stationary method (e.g. Richardson's iteration) is capable
of being constructed from d phi/dt = coefficient (laplacian phi) to
solve laplacian phi = 0. A fictitious time step of size dt = 1 along
with
a forward euler time advancement scheme would need to be used.

dpe

[Molemaker, Maarten J.]

The equation solved currently in the code is:

laplace(phi) = div(u*). The right-hand side to the laplace equation is
the divergence of the intermediate flow field. This is part of the projection
onto a divergence free velocity field. You can off course change that to have
the poisson solver solve for phi for a zero right-hand side. However, with the
current boundary conditions (homogenous neumann), the result would be fairly boring.... i.e. zero.

[Jonathan Cohen]

There are two solvers for the Poisson equation:
\nabla^2 \phi = rhs
where you are free to set rhs however you want.  As Jeroen pointed out, often you set the rhs to be the divergence of the velocity field, which results in a projection method for divergence-free flow.  But you can solve any equation, including a Laplace equation with a zero rhs.  Of course, you need interesting boundary conditions or you will just get the trivial solution (zero).

One uses a multigrid method:
http://code.google.com/p/opencurrent/source/browse/src/ocuequation/sol_mgpressure3d.h
This solver is restricted to boundary conditions at the external edges of the grid only, but is very fast.

The other uses a matrix-free preconditioned conjugate gradient (PCG) method:
http://code.google.com/p/opencurrent/source/browse/src/ocuequation/sol_pcgpressure3d.h
see the 'pcg' app for an example of how to use it:
http://code.google.com/p/opencurrent/source/browse/src/apps/pcg.cpp

The PCG solver currently solves variable-coefficitient Poisson equations:
div( A grad (\phi)) = rhs
which is a bit more general (if A is uniformly 1, this becomes just the equation above, \nabla^2 \phi = rhs).
It could be extended to handle complex boundary conditions in a pretty straightforward way, but this is not currently implemented.  If you're interested in doing this, I'd be happy to help.  This would be generally very useful in a lot of different contexts, including stationary airflow problems like you mentioned.

The PCG solver isn't in the 1.0 release, so to get that code you'll have to check the project out using hg
(instructions here: http://code.google.com/p/opencurrent/source/checkout)

-Jon

[tipara]

If I understand correctly, in the PCG solver, A is the grid called
"_coefficient" in sol_pcgpressure3d ? Also, I was wondering : is the
convergence computer-dependant ? Because I can't get the pcg app to
converge under 10^-4 within 100 iterations on mine, is it normal ?

For my problem, I need the the velocity potential on the sides of an
obstacle in a stationary flow, with Neumann conditions for x,y,z neg/
pos, and forcing dphi/dn = 0 on the sides of the obstacle. In
sol_pcgpressure3d.cu, I think the diag_preconditioner and
apply_laplacian functions should be modified to handle the "flags"
grid (1 if the obstacle is there, 0 if not), maybe simply using the
modified laplacian from the first reply ? But there must be other
functions to change that I haven't understood, mainly do_pcg, right ?

[Molemaker, Maarten J.]

A lack of convergence may be due to the insolvability of your system. For homogenous neumann boundary conditions (dphi/dn = 0), the integral of the right hand side must be discretely zero. If your solver initially appears to converge, but then seems 'stuck' at, in your case 1e-4, this would be a likely explanation. If the solver is converging monotonically, but very slowly, the condition number of your preconditioned system is bad. In most cases, this would point to a poor preconditioning.

-----Original Message-----
From: opencurr...@googlegroups.com on behalf of tipara
Sent: Sun 2/14/2010 4:40 PM
To: opencurrent-users
Subject: [opencurrent-users] Re: Working between two boundaries