Hi Jose,
Following my previous email, just some update, the epsilon actually does not
change in my shell matrix to compute the finite difference Jacobian. I don't
have an iterative procedure in the shell matrix. It simply computes the
non-linear residuals twice with perturbed flow variables. I have attached my
code. Any suggestions are welcome.
Thanks,
Feng
________________________________
From: Jose E. Roman <[email protected]>
Sent: 25 August 2022 9:45
To: feng wang <[email protected]>
Cc: petsc-users <[email protected]>
Subject: Re: [petsc-users] Slepc Question, shell-matrix
Probably the explanation is the following. Your shell matrix runs an iterative
procedure that computes the action of the matrix to a certain accuracy. This is
equivalent to multiplying by the exact matrix plus a perturbation, A+E. But
every time you apply the shell matrix, the perturbation E is slightly
different, so the eigensolver does not really see a constant matrix, it changes
at every iteration.
This also happens in nested iterations, such as inexact shift-and-invert, where
e.g. GMRES is run within each iteration of the eigensolver. For this to work,
the GMRES tolerance should be smaller than the tolerance of the eigensolver.
Similarly, if your shell matrix runs an iterative procedure, you have to use a
tolerance that is more stringent than the one used in the eigensolver.
Hope this helps.
Jose
> El 25 ago 2022, a las 11:34, feng wang <[email protected]> escribió:
>
> Hi Jose,
>
> Thanks for your reply. I have fixed one thing in my shell matrix. Now the
> conversion seems working.
>
> If I use the converted dense matrix and give it to EPSSetOperators, it
> actually converges!
>
> 50 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (7.73051721e-08)
> 51 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (5.42681175e-08)
> 52 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (2.48352708e-08)
> 53 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (1.76912430e-08)
> 54 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (2.07480734e-08)
> 55 EPS nconv=0 first unconverged value (error) 2.89521e-06+2.13473e-06i
> (1.13588981e-08)
> 56 EPS nconv=2 first unconverged value (error) 1.10143e-06+9.80005e-06i
> (8.32889697e-06)
> Solution method: krylovschur
>
> Number of requested eigenvalues: 1
> Linear eigensolve converged (2 eigenpairs) due to CONVERGED_TOL; iterations
> 56
> ---------------------- --------------------
> k ||Ax-kx||/||kx||
> ---------------------- --------------------
> 0.000003+0.000002i 7.90447e-09
> 0.000003-0.000002i 7.90447e-09
> ---------------------- --------------------
>
> If I use the shell matrix directly, the relative norm is still large. This
> really puzzles me......
>
> Thanks,
> Feng
>
> From: Jose E. Roman <[email protected]>
> Sent: 25 August 2022 7:21
> To: feng wang <[email protected]>
> Cc: petsc-users <[email protected]>
> Subject: Re: [petsc-users] Slepc Question, shell-matrix
>
> This works for simple examples such as ex9.c
> It may be an indication that there is a problem with your shell matrix.
>
> Jose
>
>
> > El 24 ago 2022, a las 19:35, feng wang <[email protected]> escribió:
> >
> > Hi Jose,
> >
> > Thanks for your reply.
> >
> > I have tried
> > PetscCall(MatConvert(slepc_A_mf,MATDENSE,MAT_INITIAL_MATRIX,&Adense)); and
> > If I do MatView(Adense,viewer), All I see was NAN. Have I missed anything
> > here?
> >
> > Thanks,
> > Feng
> >
> > From: Jose E. Roman <[email protected]>
> > Sent: 24 August 2022 16:23
> > To: feng wang <[email protected]>
> > Cc: petsc-users <[email protected]>
> > Subject: Re: [petsc-users] Slepc Question, shell-matrix
> >
> >
> >
> > > El 24 ago 2022, a las 17:03, feng wang <[email protected]> escribió:
> > >
> > > Hi Jose,
> > >
> > > Thanks for your reply.
> > >
> > > I have tried -eps_view_mat0 binary:amatrix.bin to save my shell matrix.
> > > It seems not saving the shell matrix and just created an empty file.
> >
> > Oops, I thought the shell matrix would be converted automatically. Try
> > converting it yourself:
> >
> > PetscCall(MatConvert(slepc_A_mf,MATDENSE,MAT_INITIAL_MATRIX,&Adense));
> >
> > Then use Adense in EPSSetOperators(). Or save it to disk and send it to me
> > so that I can give it a try.
> >
> > >
> > > Besides, I am trying to understand the output of slepc. I can set
> > > -eps_tol to a low value, but the final relative residual norm is still
> > > high , what does this tell me?
> >
> > I don't know. It should not happen.
> > Jose
> >
> > >
> > > Best regards,
> > > Feng
> > >
> > > From: Jose E. Roman <[email protected]>
> > > Sent: 23 August 2022 11:06
> > > To: feng wang <[email protected]>
> > > Cc: petsc-users <[email protected]>
> > > Subject: Re: [petsc-users] Slepc Question, shell-matrix
> > >
> > > You can try the following. Save your shell matrix to disk with the option
> > > -eps_view_mat0 binary:amatrix.bin
> > > then repeat the computation with ex4.c loading this file.
> > >
> > > Note that this should be done for a small problem size, because the
> > > conversion from a shell matrix implies a matrix-vector product per each
> > > column.
> > >
> > > Jose
> > >
> > > > El 23 ago 2022, a las 12:57, feng wang <[email protected]> escribió:
> > > >
> > > > Hi Jose,
> > > >
> > > > Thanks for your reply.
> > > >
> > > > It represents a linear operator. In my shell matrix, I am computing the
> > > > non-linear residuals twice with perturbed flow variables. The
> > > > matrix-vector product is computed as:
> > > >
> > > > A*v = (R(q+eps*v) - R(q-eps*v))/(2*eps)
> > > >
> > > > R is the non-linear residual. q is my flow variable and it does not
> > > > change. eps is the perturbation. A is my Jacobian matrix. Besides, for
> > > > some background, I am computing a steady RANS flow with finite volume
> > > > method and trying to do a global stability analysis by looking at the
> > > > Jacobian matrix.
> > > >
> > > > Thanks,
> > > > Feng
> > > >
> > > > From: Jose E. Roman <[email protected]>
> > > > Sent: 23 August 2022 10:32
> > > > To: feng wang <[email protected]>
> > > > Cc: petsc-users <[email protected]>
> > > > Subject: Re: [petsc-users] Slepc Question, shell-matrix
> > > >
> > > > The relative residual norms that are printed at the end are too large.
> > > > For NHEP problems, they should be below the tolerance. Don't know what
> > > > is happening. Does your shell matrix represent a linear (constant)
> > > > operator? Or does it change slightly depending on the input vector?
> > > >
> > > > > El 23 ago 2022, a las 12:14, feng wang <[email protected]>
> > > > > escribió:
> > > > >
> > > > > Hi Jose,
> > > > >
> > > > > I think the previous problem comes from my side. I have some
> > > > > uninitialized values in my part of code to compute the non-linear
> > > > > residuals. so, it produces a NAN when it tries to compute the
> > > > > matrix-vector product using finite difference. This might make the
> > > > > slepc/pestc do unexpected things.
> > > > >
> > > > > Now It seems I've got slepc running. eps_nev is set to 3 and I am
> > > > > trying to compute the ones with the largest amplitudes. Below is the
> > > > > slepc output.
> > > > >
> > > > > 14 EPS converged value (error) #0 -0.000164282 (5.36813206e-09)
> > > > > 16 EPS converged value (error) #1 -0.000160691+2.17113e-05i
> > > > > (3.37429620e-09)
> > > > > 16 EPS converged value (error) #2 -0.000160691-2.17113e-05i
> > > > > (3.37429620e-09)
> > > > > Solution method: krylovschur
> > > > >
> > > > > Number of requested eigenvalues: 2
> > > > > Linear eigensolve converged (3 eigenpairs) due to CONVERGED_TOL;
> > > > > iterations 16
> > > > > ---------------------- --------------------
> > > > > k ||Ax-kx||/||kx||
> > > > > ---------------------- --------------------
> > > > > -0.000164 0.0613788
> > > > > -0.000161+0.000022i 0.0773339
> > > > > -0.000161-0.000022i 0.0774536
> > > > > ---------------------- --------------------
> > > > >
> > > > > The values in the brackets are the absolute error (I believe) and
> > > > > they seem very low. The relative error seems quite large. Could you
> > > > > please comment on this?
> > > > >
> > > > >
> > > > > Best regards,
> > > > > Feng
> > > > >
> > > > > From: Jose E. Roman <[email protected]>
> > > > > Sent: 23 August 2022 5:24
> > > > > To: feng wang <[email protected]>
> > > > > Cc: petsc-users <[email protected]>
> > > > > Subject: Re: [petsc-users] Slepc Question, shell-matrix
> > > > >
> > > > > Please always respond to the list, otherwise the thread appears as
> > > > > unresolved in the archives of the mailing list.
> > > > >
> > > > >
> > > > > > El 22 ago 2022, a las 22:45, feng wang <[email protected]>
> > > > > > escribió:
> > > > > >
> > > > > > Hi Jose,
> > > > > >
> > > > > > I think I might have solved my problem. I have some uninitialized
> > > > > > values in my part of code to compute the right hand side. so it
> > > > > > produces a NAN when it tries to compute the matrix-vector product.
> > > > > >
> > > > > > Many thanks for your help!
> > > > > >
> > > > > > Best regards,
> > > > > > Feng
> > > > > > From: Jose E. Roman <[email protected]>
> > > > > > Sent: 22 August 2022 19:32
> > > > > > To: feng wang <[email protected]>
> > > > > > Cc: [email protected] <[email protected]>
> > > > > > Subject: Re: [petsc-users] Slepc Question, shell-matrix
> > > > > >
> > > > > > This is very strange. This error appears when the solver employs
> > > > > > the B-inner product, but in your case you don't have a B matrix, so
> > > > > > you should never see that error. Try running under valgrind to see
> > > > > > if it gives more hints.
> > > > > >
> > > > > > Jose
> > > > > >
> > > > > >
> > > > > > > El 22 ago 2022, a las 20:45, feng wang <[email protected]>
> > > > > > > escribió:
> > > > > > >
> > > > > > > Hello,
> > > > > > >
> > > > > > > I am new to Slepc and trying to work out the eigenvalues and
> > > > > > > eigenvectors of my Jacobian matrix. I am using a shell matrix to
> > > > > > > work out the matrix-vector product and I am using the default
> > > > > > > Krylov-schur method.
> > > > > > >
> > > > > > > My first attempt was not successful and I got the following
> > > > > > > errors:
> > > > > > >
> > > > > > > [0]PETSC ERROR: --------------------- Error Message
> > > > > > > --------------------------------------------------------------
> > > > > > > [0]PETSC ERROR: Missing or incorrect user input
> > > > > > > [0]PETSC ERROR: The inner product is not well defined: indefinite
> > > > > > > matrix
> > > > > > > [0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble
> > > > > > > shooting.
> > > > > > > [0]PETSC ERROR: Petsc Release Version 3.17.4, unknown
> > > > > > > [0]PETSC ERROR: cfdtest on a arch-debug named ming by feng Mon
> > > > > > > Aug 22 19:21:41 2022
> > > > > > > [0]PETSC ERROR: Configure options --with-cc=mpicc
> > > > > > > --with-cxx=mpicxx --with-fc=0 PETSC_ARCH=arch-debug
> > > > > > > [0]PETSC ERROR: #1 BV_SafeSqrt() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/include/slepc/private/bvimpl.h:130
> > > > > > > [0]PETSC ERROR: #2 BV_SquareRoot_Default() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/include/slepc/private/bvimpl.h:365
> > > > > > > [0]PETSC ERROR: #3 BVOrthogonalizeCGS1() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:101
> > > > > > > [0]PETSC ERROR: #4 BVOrthogonalizeGS() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:177
> > > > > > > [0]PETSC ERROR: #5 BVOrthonormalizeColumn() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvorthog.c:402
> > > > > > > [0]PETSC ERROR: #6 BVMatArnoldi() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/sys/classes/bv/interface/bvkrylov.c:91
> > > > > > > [0]PETSC ERROR: #7 EPSSolve_KrylovSchur_Default() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/eps/impls/krylov/krylovschur/krylovschur.c:261
> > > > > > > [0]PETSC ERROR: #8 EPSSolve() at
> > > > > > > /home/feng/cfd/slepc-3.17.1/src/eps/interface/epssolve.c:147
> > > > > > > [0]PETSC ERROR: #9 slepc_eigen_comp() at
> > > > > > > domain/cfd/slepc_eigen_solve.cpp:77
> > > > > > >
> > > > > > > Could someone please shine some light on this? I have also
> > > > > > > attached my code. The code is part of my big library and I cannot
> > > > > > > attach the whole code, sorry about this. but I am happy to
> > > > > > > provide more information. The attached code has some arrangements
> > > > > > > for halo exchange, but for the moment, it assumes it is a serial
> > > > > > > run.
> > > > > > >
> > > > > > > Many thanks,
> > > > > > > Feng
> > > > > > > <slepc_eigen_solve.cpp>
using namespace std;
# include <iomanip> // std::setprecision
# include <domain/cfd/domain.h>
#ifdef PETSC
PetscErrorCode cFdDomain::slepc_eigen_comp()
{
cout << "slepc_eigen_solve\n";
PetscErrorCode ierr;
Mat slepc_A_mf; /* eigenvalue problem matrix */
Int iq, iqs, iqe;
PetscInt blocksize;
MPI_Comm* A_COMM_WORLD;
PetscInt maxneig=8;
PetscInt nev;
EPS eps; /* eigenproblem solver context */
EPSType type;
PetscViewer viewer;
movegrid();
frame();
weights();
cout << "show vtk file\n";
vtk();
bcs();
// bcs_z(0);
dof->range( dev->getrank(), &iqs,&iqe );
////compute mean flow rhs
grad( q,dqdx );
gradb( q,dqdx );
setv( 0,dof->size(), nv, 0., rhs );
gtrhf( rhs );
gtrhm( rhs );
gtrhv( rhs );
A_COMM_WORLD = dev->getcomm();
dof->range( dev->getrank(), &iqs,&iqe );
//number of ghost cells and their global indexes (in the new numbering)
ig_mat = dof->get_igmat();
ighost = new Int [dof->size()];
ighost_local = new Int [dof->size()];
ilocal = new Int [dof->size()];
nghost=0;
nlocal=0;
for(iq=0; iq<dof->size(); iq++)
{
if(iq>=iqs && iq<iqe)
{
//non-halo cells
ilocal[nlocal] = iq;
nlocal++;
}
else
{
ighost[nghost] = ig_mat[iq];
ighost_local[nghost] = iq;
nghost++;
}
}
blocksize = nv;
//create shell matrix
PetscCall(MatCreateShell(*A_COMM_WORLD, nlocal*blocksize, nlocal*blocksize, PETSC_DETERMINE, PETSC_DETERMINE, this, &slepc_A_mf));
PetscCall(MatShellSetOperation(slepc_A_mf,MATOP_MULT,(void(*)(void))mymult_slepc));
// ierr = MatCreateBAIJ(*A_COMM_WORLD, blocksize, nlocal*blocksize, nlocal*blocksize, PETSC_DETERMINE, PETSC_DETERMINE,
// maxneig, NULL, maxneig, NULL, &petsc_A_pre); CHKERRQ(ierr);
// ierr = MatSetOption(petsc_A_pre, MAT_STRUCTURALLY_SYMMETRIC, PETSC_TRUE); CHKERRQ(ierr);
// Mat Adense;
//cout << "convert shell to dense matrix\n";
// PetscCall(MatConvert(slepc_A_mf,MATDENSE,MAT_INITIAL_MATRIX,&Adense));
//cout << "convert shell to dense matrix, done...\n";
// Create eigensolver context
PetscCall(EPSCreate(*A_COMM_WORLD,&eps));
// Set operators. In this case, it is a standard eigenvalue problem
PetscCall(EPSSetOperators(eps,slepc_A_mf,PETSC_NULL));
//PetscCall(EPSSetOperators(eps,Adense,PETSC_NULL)); //converges, works
PetscCall(EPSSetProblemType(eps,EPS_NHEP));
// Ask for the rightmost eigenvalues
PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
// Set other solver options at runtime
PetscCall(EPSSetFromOptions(eps));
// Solve the eigensystem
PetscCall(EPSSolve(eps));
/*
Optional: Get some information from the solver and display it
*/
PetscCall(EPSGetType(eps,&type));
PetscCall(PetscPrintf(*A_COMM_WORLD," Solution method: %s\n\n",type));
PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
PetscCall(PetscPrintf(*A_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* show detailed info unless -terse option is given by user */
PetscCall(PetscViewerASCIIGetStdout(*A_COMM_WORLD,&viewer));
PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
PetscCall(EPSConvergedReasonView(eps,viewer));
PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
PetscCall(PetscViewerPopFormat(viewer));
PetscCall(EPSDestroy(&eps));
PetscCall(MatDestroy(&slepc_A_mf));
ierr=0;
return ierr;
}
PetscErrorCode cFdDomain::mymult_slepc(Mat m ,Vec x, Vec y)
{
PetscErrorCode ierr;
void *ctx;
cFdDomain *myctx;
Vec localin;
MatShellGetContext(m, &ctx);
myctx = (cFdDomain*)ctx;
// PetscCall(VecGhostUpdateBegin(x, INSERT_VALUES, SCATTER_FORWARD));
// PetscCall(VecGhostUpdateEnd(x, INSERT_VALUES, SCATTER_FORWARD));
// PetscCall(VecGhostGetLocalForm(x,&localin));
ierr = myctx->petsc_set_fd_dq(x); CHKERRQ(ierr);
ierr = myctx->slepc_fd_rhs2(y); CHKERRQ(ierr);
//ierr = myctx->slepc_lin_flux(y); CHKERRQ(ierr);
//ierr = myctx->slepc_jac(x,y); CHKERRQ(ierr);
// VecGhostRestoreLocalForm(x,&localin);
return ierr;
}
PetscErrorCode cFdDomain::slepc_fd_rhs2(Vec var)
{
Int iqs, iqe;
Int iv, iq, jq;
Real delta_rhs;
PetscErrorCode ierr;
PetscScalar *array;
Real eps;
Real dq_norm,csv_norm;
dof->range( dev->getrank(), &iqs,&iqe );
fld->auxv( iqs,iqe, q,aux );
fld->cnsv( iqs,iqe, q, aux, csv);
dq_norm = 0;
csv_norm = 0;
for(iv=0; iv<nv; iv++)
{
for(iq=iqs; iq<iqe; iq++)
{
dq_norm += dq[iv][iq]*dq[iv][iq];
csv_norm += csv[iv][iq]*csv[iv][iq];
}
}
eps = sqrt(1+sqrt(csv_norm))*1.e-3/sqrt(dq_norm);
//cout << eps << " " << sqrt(csv_norm) << " " << dq_norm << "\n";
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<nq; iq++)
{
daux[iv][iq] = q[iv][iq];
}
}
//compute rhs with -eps
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<nq; iq++)
{
csv[iv][iq] = csv[iv][iq] - dq[iv][iq]*eps;
}
}
fld->csv_to_pri(iqs,iqe,csv,q);
grad( q,dqdx );
gradb( q,dqdx );
setv( 0,dof->size(), nv, 0., rhs );
gtrhf( rhs );
gtrhm( rhs );
gtrhv( rhs );
///reset q
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<dof->size(); iq++)
{
q[iv][iq] = daux[iv][iq];
}
}
//compute rhs with +eps
fld->auxv( iqs,iqe, q,aux );
fld->cnsv( iqs,iqe, q, aux, csv);
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<nq; iq++)
{
csv[iv][iq] = csv[iv][iq] + dq[iv][iq]*eps;
}
}
fld->csv_to_pri(iqs,iqe,csv,q);
grad( q,dqdx );
gradb( q,dqdx );
setv( 0,dof->size(), nv, 0., wrkq );
gtrhf( wrkq );
gtrhm( wrkq );
gtrhv( wrkq );
//matrix-vector product
ierr = VecGetArray(var,&array); CHKERRQ(ierr);
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<nlocal; iq++)
{
jq = ilocal[iq];
array[iv + nv*iq] = (wrkq[iv][jq] - rhs[iv][jq])/(2*eps);
// cout << array[iv + nv*iq] << " " << iv << " " << jq << " " << wrkq[iv][jq] << " " << rhs[iv][jq] << "\n";
}
}
ierr = VecRestoreArray(var,&array); CHKERRQ(ierr);
//reset q
for(iv=0; iv<nv; iv++)
{
for(iq=0; iq<dof->size(); iq++)
{
q[iv][iq] = daux[iv][iq];
}
}
ierr = 0;
return ierr;
}
#endif
