[petsc-dev] problem when creating a new ksp that calls another one

Wed, 15 Jul 2015 10:05:51 -0700 

Hello,

I want to make a new solver inside petsc.
This solver is based on a 2 stage method. So I have one outer iteration and one inner iteration. For the inner iteration I can use gmres or one of its variants. I succeeded to make my algorithm by creating a new ksp (that copies the linear system). So I have tried to do that without creating a new ksp. It works well with ksp examples that don't use dmda and with snes examples. 


So I have tried to make a really simple ksp (only to test) in which I 
simply call another ksp.

Here is the part of the code (and attached is the complete code):
   PetscErrorCode ierr;
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   KSPSetType(ksp,KSPFGMRES);
   ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
   PetscFunctionReturn(0);


For example with the ksp example 15:

 time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 3    ./ex15 -m 200 
-n 200 -ksp_type test  -ksp_rtol 1e-10 -pc_type mg -ksp_monitor

there is no problem, I obtain the good result.

But with the ksp example 45 (based on dmda):

time $PETSC_DIR/arch-linux2-c-debug/bin/mpirun  -np 1    ./ex45 
-ksp_type test  -ksp_rtol 1e-10 -pc_type mg -ksp_monitor


It crashes with that (I didn't change anything in ex45)
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR:  Vec is locked read only, argument # 1

[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html 
for trouble shooting.

[0]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015

[0]PETSC ERROR: ./ex45 on a arch-linux2-c-debug named extinction by 
couturie Wed Jul 15 18:45:39 2015
[0]PETSC ERROR: Configure options --download-openmpi --download-hypre 
--download-f2cblaslapack --with-fc=0
[0]PETSC ERROR: #1 VecGetArray() line 1646 in 
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #2 VecGetArray3d() line 2411 in 
/home/couturie/work/petsc-3.6.0/src/vec/vec/interface/rvector.c
[0]PETSC ERROR: #3 DMDAVecGetArray() line 77 in 
/home/couturie/work/petsc-3.6.0/src/dm/impls/da/dagetarray.c
[0]PETSC ERROR: #4 ComputeRHS() line 84 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c
[0]PETSC ERROR: #5 KSPSetUp() line 277 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #6 KSPSolve() line 546 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #7 KSPSolve_TEST() line 43 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/impls/gmres/test.c
[0]PETSC ERROR: #8 KSPSolve() line 604 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/interface/itfunc.c
[0]PETSC ERROR: #9 main() line 51 in 
/home/couturie/work/petsc-3.6.0/src/ksp/ksp/examples/tutorials/ex45.c




I don't understand since many examples in snes are good with my code and 
they use dmda.

Do any of you have an idea of the problem?


Of course, I have a bad solution that consists in creating a new ksp but 
I would like to do without, if possible.


Thank you in advance,
Regards,

Raphaël Couturier

#include <petsc/private/kspimpl.h>        /*I "petscksp.h" I*/



#undef __FUNCT__
#define __FUNCT__ "KSPSetUp_TEST"
static PetscErrorCode KSPSetUp_TEST(KSP ksp)
{
  PetscErrorCode ierr;
  PetscBool      diagonalscale;

  PetscFunctionBegin;
  ierr = PCGetDiagonalScale(ksp->pc,&diagonalscale);CHKERRQ(ierr);
  if (diagonalscale) SETERRQ1(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"Krylov method %s does not support diagonal scaling",((PetscObject)ksp)->type_name);
  ierr = KSPSetWorkVecs(ksp,9);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}












//good
#undef __FUNCT__
#define __FUNCT__ "KSPSolve_TEST"
 PetscErrorCode KSPSolve_TEST(KSP ksp) {


   PetscErrorCode ierr;
   ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
   
   KSPSetType(ksp,KSPFGMRES);
   
   
   ierr = KSPSolve(ksp, NULL, NULL); CHKERRQ(ierr);
   
   PetscFunctionReturn(0);
   








}












#undef __FUNCT__
#define __FUNCT__ "KSPCreate_TEST"
PETSC_EXTERN PetscErrorCode KSPCreate_TEST(KSP ksp)
{
  PetscErrorCode ierr;

  PetscFunctionBegin;
  ksp->data = (void*)0;

  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_PRECONDITIONED,PC_LEFT,2);CHKERRQ(ierr);
  ierr = KSPSetSupportedNorm(ksp,KSP_NORM_UNPRECONDITIONED,PC_RIGHT,1);CHKERRQ(ierr);

  ksp->ops->setup          = KSPSetUp_TEST;
  ksp->ops->solve          = KSPSolve_TEST;
  ksp->ops->destroy        = KSPDestroyDefault;
  ksp->ops->buildsolution  = KSPBuildSolutionDefault;
  ksp->ops->buildresidual  = KSPBuildResidualDefault;
  ksp->ops->setfromoptions = 0;
  ksp->ops->view           = 0;
#if defined(PETSC_USE_COMPLEX)
  SETERRQ(PetscObjectComm((PetscObject)ksp),PETSC_ERR_SUP,"This is not supported for complex numbers");
#endif
  PetscFunctionReturn(0);
}






































					Reply via email to