<div dir="ltr"><div>Hi!<br><br>Just a second try in case the first email was unintentionally missed. In short, I'm wondering if it is possible to use TSSolve without using any DM structure:<br><br>When I was debugging my program I tried to analyse the difference between the TSSolve-based time iteration program, and a version where I do the time-stepping manually using a for-loop which repeateed KSPSolve functions. The debugger showed that in the TSSolve program a pointer named ts->snes->ksp->pc->dm was nonzero, but in the for-loop program it was zero. In the TSSolve-program, and pc->dm != 0 made it enter a block of code within PCSetUp_MG that called DMCoarsen, which then crashed. So I guess I'm wondering how to tell PETSc that I'm not using any DM structure.<br>
<br></div>Best regards<br>Torquil Sørensen<br></div><div class="gmail_extra"><br><br><div class="gmail_quote">On 14 December 2013 14:46, Torquil Macdonald Sørensen <span dir="ltr"><<a href="mailto:torquil@gmail.com" target="_blank">torquil@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi!<br>
<br>
I'm getting segfault problems when trying to use TSStep or TSSolve with<br>
TSBEULER and TSLINEAR, together with PCMG. So I have reduced my code to<br>
a minimal (and mathematically trivial) test case which displays the same<br>
error, and I'm hoping someone here will spot any errors that I've made<br>
in my setup of TS.<br>
<br>
I can use my PCMG restriction and interpolation matrices with success in<br>
a single linear solve with KSPSolve (and can see that MG is being used<br>
by adding -ksp_view), but for some reason I am not able to get PCMG<br>
working with TSStep. So I'm wondering if I'm forgetting a step when when<br>
I'm setting up the TS?<br>
<br>
My test case tries to solve Id*dU/dt = Id*U, where Id is a 5x5 identity<br>
matrix, and U is a 5-vector unknown, with initial condition U = 0. The<br>
PCMG construction uses an additional coarse grid of 3 points. The fine<br>
grid vertices are labelled 0, 1, 2, 3, 4, while the coarse grid is<br>
labelled 0, 1, 2, The coarse grid vertices coincide with the fine grid<br>
vertices 0, 2 and 4. The resulting interpolation and restriction<br>
matrices arise as interpolations matrices between the FEM spaces defined<br>
using piecewise linear shape functions on these grids.<br>
<br>
The GDB message I'm getting is:<br>
<br>
Program received signal SIGSEGV, Segmentation fault.<br>
0x0000000000000000 in ?? ()<br>
(gdb) where<br>
#0 0x0000000000000000 in ?? ()<br>
#1 0x00007fa0e6bbd172 in DMCoarsen (dm=0x28e9170, comm=0x7fa0e638d9c0<br>
<ompi_mpi_comm_null>, dmc=0x295a490) at dm.c:1811<br>
#2 0x00007fa0e6ff1c1e in PCSetUp_MG (pc=0x28d71a0) at mg.c:593<br>
#3 0x00007fa0e72eb6ce in PCSetUp (pc=0x28d71a0) at precon.c:890<br>
#4 0x00007fa0e6e6f7e8 in KSPSetUp (ksp=0x28e1a50) at itfunc.c:278<br>
#5 0x00007fa0e6e70a30 in KSPSolve (ksp=0x28e1a50, b=0x297bea0,<br>
x=0x2981970) at itfunc.c:399<br>
#6 0x00007fa0e6e8ae0c in SNESSolve_KSPONLY (snes=0x28c7710) at ksponly.c:44<br>
#7 0x00007fa0e757b642 in SNESSolve (snes=0x28c7710, b=0x0, x=0x296c7d0)<br>
at snes.c:3636<br>
#8 0x00007fa0e76280cc in TSStep_Theta (ts=0x2894220) at theta.c:183<br>
#9 0x00007fa0e766a202 in TSStep (ts=0x2894220) at ts.c:2507<br>
#10 0x0000000000403363 in main (argc=3, argv=0x7fff08015378) at<br>
/home/tmac/research/schrodinger/source/testcase.cpp:97<br>
<br>
My testcase code is:<br>
<br>
#include "petscts.h"<br>
<br>
int main(int argc, char ** argv)<br>
{<br>
PetscErrorCode ierr = PetscInitialize(&argc, &argv, 0, 0);<br>
CHKERRQ(ierr);<br>
<br>
int const nb_dof_fine = 5;<br>
int const nb_dof_coarse = 3;<br>
<br>
Vec x; // A vector to hold the solution on the fine grid<br>
ierr = VecCreateSeq(PETSC_COMM_SELF, nb_dof_fine, &x); CHKERRQ(ierr);<br>
<br>
Mat Id; // Identity matrix on the vector space of the fine grid<br>
ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nb_dof_fine, nb_dof_fine, 1,<br>
0, &Id); CHKERRQ(ierr);<br>
for(int i = 0; i != nb_dof_fine; i++) { MatSetValue(Id, i, i, 1.0,<br>
INSERT_VALUES); }<br>
ierr = MatAssemblyBegin(Id, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(Id, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatView(Id, PETSC_VIEWER_STDOUT_SELF);<br>
<br>
Mat interp; // Multigrid interpolation matrix<br>
ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nb_dof_fine, nb_dof_coarse,<br>
2, 0, &interp); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 0, 0, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 1, 0, 0.5, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 1, 1, 0.5, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 2, 1, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 3, 1, 0.5, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 3, 2, 0.5, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(interp, 4, 2, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatAssemblyBegin(interp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(interp, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatView(interp, PETSC_VIEWER_STDOUT_SELF);<br>
<br>
Mat restrict; // Multigrid restriction matrix, not the transpose of<br>
"interp"<br>
ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nb_dof_coarse, nb_dof_fine,<br>
1, 0, &restrict); CHKERRQ(ierr);<br>
ierr = MatSetValue(restrict, 0, 0, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(restrict, 1, 2, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatSetValue(restrict, 2, 4, 1.0, INSERT_VALUES); CHKERRQ(ierr);<br>
ierr = MatAssemblyBegin(restrict, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatAssemblyEnd(restrict, MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);<br>
ierr = MatView(restrict, PETSC_VIEWER_STDOUT_SELF);<br>
<br>
{ // Solving the linear problem Id*x = 0 works<br>
KSP ksp;<br>
ierr = KSPCreate(PETSC_COMM_SELF, &ksp); CHKERRQ(ierr);<br>
ierr = KSPSetType(ksp, KSPFGMRES); CHKERRQ(ierr);<br>
ierr = KSPSetOperators(ksp, Id, Id, SAME_PRECONDITIONER);<br>
CHKERRQ(ierr);<br>
<br>
PC pc;<br>
ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);<br>
ierr = PCSetType(pc, PCMG); CHKERRQ(ierr);<br>
ierr = PCMGSetLevels(pc, 2, 0); CHKERRQ(ierr);<br>
ierr = PCMGSetGalerkin(pc, PETSC_TRUE); CHKERRQ(ierr);<br>
ierr = PCMGSetInterpolation(pc, 1, interp); CHKERRQ(ierr);<br>
ierr = PCMGSetRestriction(pc, 1, restrict); CHKERRQ(ierr);<br>
<br>
ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);<br>
<br>
Vec b; // Trivial RHS for the linear problem<br>
ierr = VecDuplicate(x, &b); CHKERRQ(ierr);<br>
ierr = VecSet(b, 0); CHKERRQ(ierr);<br>
ierr = KSPSolve(ksp, b, x); CHKERRQ(ierr);<br>
ierr = VecView(x, PETSC_VIEWER_STDOUT_SELF);<br>
<br>
ierr = VecDestroy(&b); CHKERRQ(ierr);<br>
ierr = KSPDestroy(&ksp); CHKERRQ(ierr);<br>
<br>
ierr = PetscPrintf(PETSC_COMM_SELF, "KSPSolve worked\n");<br>
CHKERRQ(ierr);<br>
}<br>
<br>
{ // This block of code causes a segfault<br>
TS ts;<br>
ierr = TSCreate(PETSC_COMM_SELF, &ts); CHKERRQ(ierr);<br>
ierr = TSSetProblemType(ts, TS_LINEAR); CHKERRQ(ierr);<br>
ierr = TSSetType(ts, TSBEULER); // Use an implicit scheme<br>
ierr = TSSetInitialTimeStep(ts, 0.0, 0.1); CHKERRQ(ierr);<br>
ierr = TSSetSolution(ts, x); CHKERRQ(ierr);<br>
ierr = TSSetRHSFunction(ts, 0, TSComputeRHSFunctionLinear, 0);<br>
CHKERRQ(ierr);<br>
ierr = TSSetRHSJacobian(ts, Id, Id,<br>
TSComputeRHSJacobianConstant, 0); CHKERRQ(ierr);<br>
ierr = TSSetIFunction(ts, 0, TSComputeIFunctionLinear, 0);<br>
CHKERRQ(ierr);<br>
ierr = TSSetIJacobian(ts, Id, Id, TSComputeIJacobianConstant, 0);<br>
<br>
KSP ksp;<br>
ierr = TSGetKSP(ts, &ksp); CHKERRQ(ierr);<br>
ierr = KSPSetType(ksp, KSPFGMRES); CHKERRQ(ierr);<br>
<br>
PC pc;<br>
ierr = KSPGetPC(ksp, &pc); CHKERRQ(ierr);<br>
ierr = PCSetType(pc, PCMG); CHKERRQ(ierr);<br>
ierr = PCMGSetLevels(pc, 2, 0); CHKERRQ(ierr);<br>
ierr = PCMGSetGalerkin(pc, PETSC_TRUE); CHKERRQ(ierr);<br>
ierr = PCMGSetInterpolation(pc, 1, interp); CHKERRQ(ierr);<br>
ierr = PCMGSetRestriction(pc, 1, restrict); CHKERRQ(ierr);<br>
<br>
ierr = TSSetUp(ts); CHKERRQ(ierr);<br>
ierr = TSStep(ts); CHKERRQ(ierr);<br>
<br>
ierr = TSDestroy(&ts); CHKERRQ(ierr);<br>
}<br>
<br>
ierr = MatDestroy(&restrict); CHKERRQ(ierr);<br>
ierr = MatDestroy(&interp); CHKERRQ(ierr);<br>
ierr = VecDestroy(&x); CHKERRQ(ierr);<br>
ierr = MatDestroy(&Id); CHKERRQ(ierr);<br>
ierr = PetscFinalize(); CHKERRQ(ierr);<br>
<br>
return 0;<br>
}<br>
<br>
The actual output of the program when running without a debugger is:<br>
<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
row 0: (0, 1)<br>
row 1: (1, 1)<br>
row 2: (2, 1)<br>
row 3: (3, 1)<br>
row 4: (4, 1)<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
row 0: (0, 1)<br>
row 1: (0, 0.5) (1, 0.5)<br>
row 2: (1, 1)<br>
row 3: (1, 0.5) (2, 0.5)<br>
row 4: (2, 1)<br>
Matrix Object: 1 MPI processes<br>
type: seqaij<br>
row 0: (0, 1)<br>
row 1: (2, 1)<br>
row 2: (4, 1)<br>
Vector Object: 1 MPI processes<br>
type: seq<br>
0<br>
0<br>
0<br>
0<br>
0<br>
KSPSolve worked<br>
[0]PETSC ERROR:<br>
------------------------------------------------------------------------<br>
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,<br>
probably memory access out of range<br>
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger<br>
[0]PETSC ERROR: or see<br>
<a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind[0]PETSC</a><br>
ERROR: or try <a href="http://valgrind.org" target="_blank">http://valgrind.org</a> on GNU/linux and Apple Mac OS X to<br>
find memory corruption errors<br>
[0]PETSC ERROR: likely location of problem given in stack below<br>
[0]PETSC ERROR: --------------------- Stack Frames<br>
------------------------------------<br>
[0]PETSC ERROR: Note: The EXACT line numbers in the stack are not available,<br>
[0]PETSC ERROR: INSTEAD the line number of the start of the function<br>
[0]PETSC ERROR: is given.<br>
[0]PETSC ERROR: [0] DMCoarsen line 1809 src/dm/interface/dm.c<br>
[0]PETSC ERROR: [0] PCSetUp_MG line 546 src/ksp/pc/impls/mg/mg.c<br>
[0]PETSC ERROR: [0] PCSetUp line 868 src/ksp/pc/interface/precon.c<br>
[0]PETSC ERROR: [0] KSPSetUp line 192 src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: [0] KSPSolve line 356 src/ksp/ksp/interface/itfunc.c<br>
[0]PETSC ERROR: [0] SNESSolve_KSPONLY line 14<br>
src/snes/impls/ksponly/ksponly.c<br>
[0]PETSC ERROR: [0] SNESSolve line 3589 src/snes/interface/snes.c<br>
[0]PETSC ERROR: [0] TSStep_Theta line 162<br>
src/ts/impls/implicit/theta/theta.c<br>
[0]PETSC ERROR: [0] TSStep line 2499 src/ts/interface/ts.c<br>
[0]PETSC ERROR: --------------------- Error Message<br>
------------------------------------<br>
[0]PETSC ERROR: Signal received!<br>
[0]PETSC ERROR:<br>
------------------------------------------------------------------------<br>
[0]PETSC ERROR: Petsc Release Version 3.4.3, Oct, 15, 2013<br>
[0]PETSC ERROR: See docs/changes/index.html for recent updates.<br>
[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.<br>
[0]PETSC ERROR: See docs/index.html for manual pages.<br>
[0]PETSC ERROR:<br>
------------------------------------------------------------------------<br>
[0]PETSC ERROR: ./testcase.elf on a linux-gnu-c-debug named lenovo by<br>
tmac Sat Dec 14 14:35:56 2013<br>
[0]PETSC ERROR: Libraries linked from /tmp/petsc-3.4.3/linux-gnu-c-debug/lib<br>
[0]PETSC ERROR: Configure run at Fri Nov 22 11:07:41 2013<br>
[0]PETSC ERROR: Configure options --with-shared-libraries<br>
--with-debugging=1 --useThreads 0 --with-clanguage=C++ --with-c-support<br>
--with-fortran-interfaces=1 --with-scalar-type=complex<br>
--with-mpi-dir=/usr/lib/openmpi --with-blas-lib=-lblas<br>
--with-lapack-lib=-llapack --with-blacs=1<br>
--with-blacs-include=/usr/include<br>
--with-blacs-lib="[/usr/lib/libblacsCinit-openmpi.so,/usr/lib/libblacs-openmpi.so]"<br>
--with-scalapack=1 --with-scalapack-include=/usr/include<br>
--with-scalapack-lib=/usr/lib/libscalapack-openmpi.so --with-mumps=1<br>
--with-mumps-include=/usr/include<br>
--with-mumps-lib="[/usr/lib/libdmumps.so,/usr/lib/libzmumps.so,/usr/lib/libsmumps.so,/usr/lib/libcmumps.so,/usr/lib/libmumps_common.so,/usr/lib/libpord.so]"<br>
--with-umfpack=1 --with-umfpack-include=/usr/include/suitesparse<br>
--with-umfpack-lib="[/usr/lib/libumfpack.so,/usr/lib/libamd.so]"<br>
--with-cholmod=1 --with-cholmod-include=/usr/include/suitesparse<br>
--with-cholmod-lib=/usr/lib/libcholmod.so --with-spooles=1<br>
--with-spooles-include=/usr/include/spooles<br>
--with-spooles-lib=/usr/lib/libspooles.so --with-ptscotch=1<br>
--with-ptscotch-include=/usr/include/scotch<br>
--with-ptscotch-lib="[/usr/lib/libptesmumps.so,/usr/lib/libptscotch.so,/usr/lib/libptscotcherr.so]"<br>
--with-fftw=1 --with-fftw-include=/usr/include<br>
--with-fftw-lib="[/usr/lib/x86_64-linux-gnu/libfftw3.so,/usr/lib/x86_64-linux-gnu/libfftw3_mpi.so]"<br>
--with-superlu=1 --with-superlu-include=/usr/include/superlu<br>
--with-superlu-lib=/usr/lib/libsuperlu.so<br>
--CXX_LINKER_FLAGS=-Wl,--no-as-needed<br>
[0]PETSC ERROR:<br>
------------------------------------------------------------------------<br>
[0]PETSC ERROR: User provided function() line 0 in unknown directory<br>
unknown file<br>
--------------------------------------------------------------------------<br>
MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD<br>
with errorcode 59.<br>
<br>
NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.<br>
You may or may not see output from other processes, depending on<br>
exactly when Open MPI kills them.<br>
--------------------------------------------------------------------------<br>
<br>
Best regards<br>
<span class="HOEnZb"><font color="#888888">Torquil Sørensen<br>
<br>
</font></span></blockquote></div><br></div>