[petsc-users] saving results

[Sal Am]

 It is a finite-element problem of an RF antenna dielectric interaction
where all the non-zero elements are on the diagonal of the sparse matrix
(if that is relevant).
More about the matrix: 25Mx25M, 2B non-zeros. I checked the KSPMonitor, it
seems that I have to write my own routine?

> Running for a Krylov method for tens of thousands of iterations is very
> rarely recommended.
>
GAMG and BCGS is the only ones that have actually worked for me so far, I
increased the max_it because it was not enough with the default one. With
the default tolerance I got ||r(i)||/||b|| in the order of 10^-3, but I
need more accuracy so I increased the tol to 10^-7 and that is when it
crashed after ~51000 iterations.

@Barry

> I'm guessing you mange your own time stepping (and nonlinear solver if
> there is one)?
>

It is the default from the solver (attached the short code).

> As Jed says it doesn't make sense to save "partial" solutions within the
> linear solver.  Just save it every several time-steps.
>

Yes maybe I worded it wrong. I did mean to save checkpoints basically such
that the next time it is running it can pick up from where it crashed.

   Also the code should not be "crashing" at seemingly long times (after
> hours) with a Segmentation fault. Send us the full error message and we'll
> see if there is some way we can help you fix this problem.
>
I have added this before, but seemingly they conclusion was that it is
memory error...although I am not sure how, but I have added the entire
error.

On Mon, Feb 18, 2019 at 8:21 PM Smith, Barry F. <bsmith at mcs.anl.gov> wrote:

>
>    I'm guessing you mange your own time stepping (and nonlinear solver if
> there is one)?
>
>    You can save the solution with a call to VecView() and then reload the
> solution with a VecLoad() but you need to manage any other restart data
> that you may need, like the value of the current time etc.
>
>    As Jed says it doesn't make sense to save "partial" solutions within
> the linear solver.  Just save it every several time-steps.
>
>    Barry
>
>    Also the code should not be "crashing" at seemingly long times (after
> hours) with a Segmentation fault. Send us the full error message and we'll
> see if there is some way we can help you fix this problem.
>
>
> > On Feb 18, 2019, at 9:47 AM, Sal Am via petsc-users <
> petsc-users at mcs.anl.gov> wrote:
> >
> > Is there a function/command line option to save the solution as it is
> solving (and read in the file from where it crashed and keep iterating from
> there perhaps)?
> > Had a seg fault and all the results until that point seems to have been
> lost.
> >
> >
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190219/980421ce/attachment-0001.html>
-------------- next part --------------
static char help[] = "Solves a complex sparse linear system in parallel with KSP.\n\n";

#include <petscksp.h>
#include <iostream>
int main(int argc,char **args){
    Vec            x,b;           /* approx solution, RHS */
    Mat            A, F;             /* linear system matrix */
    KSP            ksp;           /* linear solver context */
    // PetscReal      norm;          /* norm of solution error */
    PC             pc;            /* petsc preconditioner object */
    PetscMPIInt    rank, size;    /* petsc mpi rank and size object */
    PetscViewer    viewer;        /* petsc viewer object needed for vis. and reading */
    PetscErrorCode ierr;

    // Initialise PETSc needed in the beginning of program, includes MPI initialisation. 
    ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
    MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
    MPI_Comm_size(PETSC_COMM_WORLD,&size);

    // Check if configured with complex
    #if !defined(PETSC_USE_COMPLEX)
      SETERRQ(PETSC_COMM_WORLD,1,"This example requires complex numbers");
    #endif
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
           Read the matrix and right-hand-side vector that defines
           the linear system, Ax = b.
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Reading vector in binary from Vector_b.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"../petscpy/petscFiles/Vector_b.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
    ierr = VecCreate(PETSC_COMM_WORLD, &b);CHKERRQ(ierr);
    ierr = VecLoad(b,viewer); CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Vector b read.\n");CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);

    ierr = PetscPrintf(PETSC_COMM_WORLD,"Reading matrix in binary from Matrix_A.dat ...\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"../petscpy/petscFiles/Matrix_A.dat",FILE_MODE_READ,&viewer);CHKERRQ(ierr);
    ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
    ierr = MatLoad(A,viewer);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Matrix A read.\n");CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    // Duplicate vector b for solution x, content not duplicated.
    ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                  Create the linear solver and set various options
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

    
    // Create linear solver context
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Creating KSPCreate\n");CHKERRQ(ierr);
    ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);

    /* Set operators. Here the matrix that defines the linear system
       also serves as the preconditioning constructor matrix.*/
    ierr = PetscPrintf(PETSC_COMM_WORLD,"KSPSetOperators initiated\n");CHKERRQ(ierr);
    ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);

    /* Set runtime options, e.g., -ksp_type <type> -pc_type <type>
     -ksp_monitor -ksp_rtol <rtol>. Here MUMPS is used if no runtime options is provided. */
    // ierr = KSPSetType(ksp, KSPGMRES);CHKERRQ(ierr);
    // ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
    // ierr = PCSetType(pc, PCLU);CHKERRQ(ierr);
    // ierr = PCFactorSetMatSolverType(pc, MATSOLVERMUMPS);CHKERRQ(ierr);
    ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);

    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                        Solve the linear system
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    /* KSPSetUp and KSPSetUpOnBlocks are optional. 
    Used to enable more precise profiling of setting up the preconditioner */ 
    ierr = PetscPrintf(PETSC_COMM_WORLD,"KSPSetUp initiated\n");CHKERRQ(ierr);
    ierr = KSPSetUp(ksp);CHKERRQ(ierr);
    // ierr = PCFactorGetMatrix(pc, &F);
    //ierr = KSPSetUpOnBlocks(ksp);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Solving the system Ax=b\n");CHKERRQ(ierr);
    ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
    ierr = PetscPrintf(PETSC_COMM_WORLD,"Done solving.\n");CHKERRQ(ierr);

    // Uncomment the following to view the solution in terminal
    // ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
    /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
                        Save solution and clean up
       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    ierr = PetscPrintf(PETSC_COMM_WORLD, "Writing solution x in PETSc binary file Vector_x_petsc.dat\n");CHKERRQ(ierr);
    ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD, "../petscpy/petscFiles/Vector_x_petsc.dat", FILE_MODE_WRITE, &viewer);CHKERRQ(ierr);
    ierr = PetscViewerBinarySetSkipHeader(viewer, PETSC_TRUE);CHKERRQ(ierr);
    ierr = VecView(x,viewer);CHKERRQ(ierr);CHKERRQ(ierr);
    ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
    /*
       Free work space.  All PETSc objects should be destroyed when they
       are no longer needed.
    */
    ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
    ierr = VecDestroy(&x);CHKERRQ(ierr);
    ierr = VecDestroy(&b);CHKERRQ(ierr); 
    ierr = MatDestroy(&A);CHKERRQ(ierr);
    ierr = PetscFinalize();
    return ierr;
}


/*TEST

   build:
      requires: complex, MUMPS
      optional: SuperLU_dist, SuperLU

   test:
      mpiexec -n 4 ./test -ksp_type gmres -pc_type lu -pc_factor_mat_solver_type mumps -ksp_max_it 1 -ksp_monitor_true_residual

TEST*/
-------------- next part --------------
A non-text attachment was scrubbed...
Name: OneStripAntenna.e42955
Type: application/octet-stream
Size: 472700 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190219/980421ce/attachment-0001.obj>