boutproject / BOUT-dev

BOUT++: Plasma fluid finite-difference simulation code in curvilinear coordinate systems
http://boutproject.github.io/
GNU Lesser General Public License v3.0
183 stars 95 forks source link

Issues with MatFDColoringSetFunction from PETSc #1160

Open ZedThree opened 6 years ago

ZedThree commented 6 years ago

With #1159, we will be clean to both -Wall and -Wextra on gcc 8.1, expect for the following two warnings (slightly abridged):

petsc.cxx: In member function ‘virtual int PetscSolver::init(int, BoutReal)’:
petsc.cxx:447:75: warning: cast between incompatible function types from \
  ‘PetscErrorCode (*)(TS, BoutReal, Vec, Vec, void*)’ ...  to \
  ‘PetscErrorCode (*)()’ ... [-Wcast-function-type]

   ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))solver_f,this);CHKERRQ(ierr);
                                                                           ^~~~~~~~

imex-bdf2.cxx: In member function ‘void IMEXBDF2::constructSNES(_p_SNES**)’:
imex-bdf2.cxx:615:69: warning: cast between incompatible function types from \
  ‘PetscErrorCode (*)(SNES, Vec, Vec, void*)’ ... to \
  ‘PetscErrorCode (*)()’ ... [-Wcast-function-type]

       MatFDColoringSetFunction(fdcoloring,(PetscErrorCode (*)(void))FormFunctionForColoring,this);

The problem is that the second argument to MatFDColoringSetFunction must be (PetscErrorCode (*)(void)), but the actual type of e.g. solver_f is PetscErrorCode (*)(TS, BoutReal, Vec, Vec, void*), hence the cast.

But these are not compatible functions in C++, hence the warnings. As far as I can tell, PETSc unconditionally casts them to PetscErrorCode (*)(void*, BoutReal, Vec, Vec, void*) anyway, so I don't know why they don't just take that instead.

The following is a horrible "fix":

ierr = MatFDColoringSetFunction(matfdcoloring, 
   reinterpret_cast<PetscErrorCode (*)(void)>(reinterpret_cast<void*>(&solver_f)),
   this);CHKERRQ(ierr);

It silences the warning by first casting to void *, which is universally compatible, hence it can be used as a buffer between what we've got and what PETSc wants. Of course, it looks horrendous and reeks to high hell.

Personally, I think PETSc should fix it their end, but even in that case, they won't backport it.


As it turns out, solver_f does have the wrong signature! It has 5 arguments instead of 4. Not sure the first parameter even has the correct type.

ZedThree commented 6 years ago

Two more things I've found:

  1. ae07134 changed MatFDColoringSetFunction to use solver_f instead of SNESTSFormFunction. Should this be reverted? Note we'd still need to do the horrible double reinterpret_cast to silence the warning.

  2. I don't think we're doing the finite difference colouring correctly, which is probably why we've not noticed passing solver_f would cause a crash if used. From the PETSc docs:

Notes: This function is usually used automatically by SNES (when one uses SNESSetJacobian() with the argument SNESComputeJacobianDefaultColor()) and only needs to be used by someone computing a matrix via coloring directly by calling MatFDColoringApply()

We never call MatFDColoringApply, or MatFDColoringSetUp.

ZedThree commented 6 years ago

This doesn't seem to be a problem because we don't hit it in most cases. To hit it, one needs to use a PETSc preconditioner. Running test-solver with -pc_type lu -J_slowfd then gives this error:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR:   
[0]PETSC ERROR: Must call TSSetRHSJacobian() and / or TSSetIJacobian()
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.9.2, May, 20, 2018 
[0]PETSC ERROR: ./test_solver on a  named fusion125 by peter Wed Jul 11 17:25:23 2018
[0]PETSC ERROR: Configure options --with-cc=mpicc --with-cxx=mpic++ --with-fc=mpif90 --with-clanguage=cxx --with-openmp=yes --with-mpi=yes --with-precision=double --with-scalar-type=real --with-shared-libraries=0 --with-debugging --with-make-np=6 --download-hypre --download-scalapack --prefix=/home/peter/Tools/petsc-3.9.2-gcc-8.1.1
[0]PETSC ERROR: #1 TSComputeRHSJacobian() line 539 in /home/peter/Tools/petsc-3.9.2/src/ts/interface/ts.c
[0]PETSC ERROR: #2 init() line 412 in petsc.cxx

I think TSSetRHSJacobian should likely be called with solver_rhsjacobian, but I'm not sure what Amat/Pmat should be -- J?