[FLASH-USERS] Multigrid boundaries with PFFT top level solver

[Klaus Weide]

Hello Joshua,

I will forward (in separate mail) a response that I wrote recently to a 
very similar question.

Some further comments below.

On Fri, 17 Jun 2016, Joshua Wall wrote:
>     I've been attempting to do a mixed boundary simulation that has
> different numbers of root grids in each direction. I think this can be
> accomplished with the multigrid solver combined with the top level PFFT
> solver (Hybrid I think its called in the user guide). But I've noticed in
> Gravity_getPotentialListOfBlocks for the MG solver that it sets all the
> boundary types equal to the type which grav_boundary_type is set to in the
> par file: [...]

Yes, but see the BHTree implementation under GravityMain/Poisson/ for an 
example of how this could be generalized to allow for specifying 
some sort of "mixed" gravity boundary conditions.

> Then later on in hgSolve it compares the boundary type thats passed with
> the type initialized from hgPFFTinit.F90 (which reads xl_mg_boundary_type,
> xr, yl, yr, etc). When it sees these are different, it throws an error 
> [...]

Yes, I believe this should happen. (At least as long as the code requires 
specifying [xyz][lr]_mg_boundary_type parameters separately from the BCs 
at the Gravity or Grid_solvePoisson levels, there should be a check for 
consistency.)

> So my question is can I safely change the setting in
> Gravity_getPotentialListOfBlocks to read the six different boundaries used
> by hgSolve and read by hgPFFTinit, since it appears it should be alright to
> solve using different ones (and its mentioned in the user guide that it can
> do mixed boundaries)? I ask because it seems almost as if someone changed
> it to only do one type of boundary because there may have been an error
> before, and I'm always a bit cautious to change such a integral piece of
> code as this.

Your approach MAY work. We have tested some combinations of mixed boundary 
conditions for the hybrid multigrid Poisson solver at the 
Grid_solvePoisson level, with unit tests that call Grid_solvePoisson 
directly (rather than via the Gravity unit).

The fact that Gravity_potentialListOfBlocks (in the implementation 
for Multigrid, at least) allows "only one type of boundary" goes back to
this being the only BCs originally implemented, and mixed BCs not being 
sufficiently tested to claim support.


> For clarity, here's my par file settings for the domain and gravity for
> what I'd like to do:
> 
>[...] 
>  39 xl_boundary_type  = "periodic" # "reflecting"
>  40 xr_boundary_type  = "periodic" # "reflecting"
>  41 yl_boundary_type  = "periodic" # "reflecting"
>  42 yr_boundary_type  = "periodic" # "reflecting"
>  43 zl_boundary_type  = "outflow" #"periodic" # "reflecting"
>  44 zr_boundary_type  = "outflow" #"periodic" # "reflecting"
> ...
> 142 grav_boundary_type = "periodic" # "isolated"
> [...]
> 146 xl_mg_boundary_type  = "periodic"
> 147 xr_mg_boundary_type  = "periodic"
> 148 yl_mg_boundary_type  = "periodic"
> 149 yr_mg_boundary_type  = "periodic"
> 150 zl_mg_boundary_type  = "outflow" #"periodic"
> 151 zr_mg_boundary_type  = "outflow" #"periodic"

My idea about how to specify what you want, instead of your line 142, goes 
like this: (with the appropriate code changes to support these extended 
runtime parameters)

  grav_boundary_type = "mixed"
  grav_boundary_type_x = "periodic"
  grav_boundary_type_y = "periodic"
  grav_boundary_type_z = "neumann" ! if you want zero-gradient potential fcn


HTH,

Klaus