[petsc-users] Problems about Assemble DMComposite Precondition Matrix
Yingjie Wu
yjwu16 at gmail.com
Mon Nov 5 08:59:31 CST 2018
Dear Petsc developer:
Hi,
I have recently studied the preconditioner of the program, and some
problems have arisen. Please help me to solve them.
At present, I have written a program to solve the system of non-linear
equations. The Matrix Free method has been used to calculate results. But I
want to add a preprocessing matrix to it.
I used the DMComposite object, which stores two sub-DM objects and a single
value (two physical field variables and one variable). I want to use
MatGetLocalSubMatrix to assign each physical field sub precondition matrix.
At the same time, because my DM object is two-dimensional, I use
MatSetValuesStencil() to fill the sub matrix.
At present, I just want to fill in a unit matrix (for global vectors) as
the precondition matrix of Matrix Free Method (just to test whether it can
be used, the unit matrix has no preprocessing effect). But the procedure
was wrong.
yjwu at yjwu-XPS-8910:~/petsc-3.10.1/src/snes/examples/tutorials$
mpiexec -n 1 ./ex216 -f wu-readtwogroups -snes_mf_operator -snes_view
-snes_converged_reason -snes_monitor -ksp_converged_reason
-ksp_monitor_true_residual
0 SNES Function norm 8.235090086536e-02
iter = 0, SNES Function norm 0.0823509
iter = 0, Keff ======= 1.
Linear solve did not converge due to DIVERGED_PCSETUP_FAILED iterations 0
PCSETUP_FAILED due to FACTOR_NUMERIC_ZEROPIVOT
Nonlinear solve did not converge due to DIVERGED_LINEAR_SOLVE iterations 0
SNES Object: 1 MPI processes
type: newtonls
maximum iterations=50, maximum function evaluations=10000
tolerances: relative=1e-08, absolute=1e-50, solution=1e-08
total number of linear solver iterations=0
total number of function evaluations=1
norm schedule ALWAYS
SNESLineSearch Object: 1 MPI processes
type: bt
interpolation: cubic
alpha=1.000000e-04
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15,
lambda=1.000000e-08
maximum iterations=40
KSP Object: 1 MPI processes
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt
Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=10000, initial guess is zero
tolerances: relative=1e-05, absolute=1e-50, divergence=10000.
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI processes
type: ilu
out-of-place factorization
0 levels of fill
tolerance for zero pivot 2.22045e-14
matrix ordering: natural
factor fill ratio given 1., needed 1.
Factored matrix follows:
Mat Object: 1 MPI processes
type: seqaij
rows=961, cols=961
package used to perform factorization: petsc
total: nonzeros=4625, allocated nonzeros=4625
total number of mallocs used during MatSetValues calls =0
not using I-node routines
linear system matrix followed by preconditioner matrix:
Mat Object: 1 MPI processes
type: mffd
rows=961, cols=961
Matrix-free approximation:
err=1.49012e-08 (relative error in function evaluation)
The compute h routine has not yet been set
Mat Object: 1 MPI processes
type: seqaij
rows=961, cols=961
total: nonzeros=4625, allocated nonzeros=4625
total number of mallocs used during MatSetValues calls =0
not using I-node routines
It seems that there are elements on the diagonal line that are not filled,
but I don't understand what went wrong. Part of the code of my program is
as follows. I added all the code to the appendix. I have tested that for
the entire Jacobian matrix assignment unit matrix (no submatrix, direct
operation of the global preprocessing matrix), the program can run
normally. However, since the problems developed later may be more complex,
involving multiple physical fields, it would be easier to implement
submatrix.
......
/* set DMComposite */
ierr =
DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da1);CHKERRQ(ierr);
ierr = DMSetFromOptions(user.da1);CHKERRQ(ierr);
ierr = DMSetUp(user.da1);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr);
ierr =
DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da2);CHKERRQ(ierr);
ierr = DMSetFromOptions(user.da2);CHKERRQ(ierr);
ierr = DMSetUp(user.da2);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr);
ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr);
......
/* in FormJacobian(SNES snes,Vec U,Mat J,Mat B,void *ctx) */
ierr =
DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr =
DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = DMCompositeGetLocalISs(user->packer,&is);CHKERRQ(ierr);
ierr = MatGetLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
ierr = MatGetLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.j = j; row.i = i;
unit = 1.0;
ierr =
MatSetValuesStencil(B11,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
}
}
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.j = j; row.i = i;
unit = 1.0;
ierr =
MatSetValuesStencil(B22,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatRestoreLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
ierr = MatRestoreLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
unit = 1.0;
row1 = 960;//last row global index
ierr = MatSetValues(B,1,&row1,1,&row1,&unit,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
......
Thanks,
Yingjie
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181105/59aceb11/attachment-0001.html>
-------------- next part --------------
static const char help[] = "Solves PDE optimization problem using full-space method, treats state and adjoint variables separately.\n\n";
/*T
requires: !single
T*/
#include <petscdm.h>
#include <petscdmda.h>
#include <petscdmredundant.h>
#include <petscdmcomposite.h>
#include <petscpf.h>
#include <petscsnes.h>
/*
w - design variables (what we change to get an optimal solution)
u - state variables (i.e. the PDE solution)
lambda - the Lagrange multipliers
U = (w u lambda)
fu, fw, flambda contain the gradient of L(w,u,lambda)
FU = (fw fu flambda)
In this example the PDE is
Uxx = 2
u(0) = w(0), thus this is the free parameter
u(1) = 0
the function we wish to minimize is
\integral u^{2}
The exact solution for u is given by u(x) = x*x - 1.25*x + .25
Use the usual centered finite differences.
Note we treat the problem as non-linear though it happens to be linear
See ex22.c for the same code, but that interlaces the u and the lambda
*/
typedef struct {
DM red1,da1,da2;
DM packer;
} UserCtx;
extern PetscErrorCode FormInitialGuess(UserCtx*,Vec,Mat,Mat);
extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
int main(int argc,char **argv)
{
PetscErrorCode ierr;
PetscInt its;
Vec U,FU,vlambda,vphi1,vphi2;
Mat A1,A2;/*Matrix for storage initial guess*/
Mat B;
SNES snes;
UserCtx user;
PetscViewer viewer;
char file[PETSC_MAX_PATH_LEN];/* input file name */
PetscBool flg;
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
/*Read the file and store them in the matrix*/
ierr = PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,&flg);CHKERRQ(ierr);
if (!flg)
{
SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f option");
}
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A1);CHKERRQ(ierr);
ierr = MatSetFromOptions(A1);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&A2);CHKERRQ(ierr);
ierr = MatSetFromOptions(A2);CHKERRQ(ierr);
ierr = MatLoad(A1,viewer);CHKERRQ(ierr);
ierr = MatLoad(A2,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
/* Create a global vector that includes a single redundant array and two da arrays */
ierr = DMCompositeCreate(PETSC_COMM_WORLD,&user.packer);CHKERRQ(ierr);
/* Set neutron flux DM */
ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da1);CHKERRQ(ierr);
ierr = DMSetFromOptions(user.da1);CHKERRQ(ierr);
ierr = DMSetUp(user.da1);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.da1);CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,20,24,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&user.da2);CHKERRQ(ierr);
ierr = DMSetFromOptions(user.da2);CHKERRQ(ierr);
ierr = DMSetUp(user.da2);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.da2);CHKERRQ(ierr);
ierr = DMRedundantCreate(PETSC_COMM_WORLD,0,1,&user.red1);CHKERRQ(ierr);
ierr = DMCompositeAddDM(user.packer,user.red1);CHKERRQ(ierr);
ierr = DMDASetFieldName(user.da1,0,"phi1");CHKERRQ(ierr);
ierr = DMDASetFieldName(user.da2,0,"phi2");CHKERRQ(ierr);
//ierr = DMDASetFieldName(user.red1,0,"lambda");CHKERRQ(ierr);
ierr = DMCreateGlobalVector(user.packer,&U);CHKERRQ(ierr);
ierr = VecDuplicate(U,&FU);CHKERRQ(ierr);
/* Form initial Guess lambda, phi1, phi2 */
ierr = FormInitialGuess(&user,U,A1,A2);CHKERRQ(ierr);
/* create nonlinear solver */
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
ierr = DMCreateMatrix(user.packer,&B);CHKERRQ(ierr);
/* This example does not correctly allocate off-diagonal blocks. These options allows new nonzeros (slow). */
ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr);
ierr = SNESSetFunction(snes,FU,FormFunction,&user);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,NULL,B,FormJacobian,&user);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = SNESMonitorSet(snes,Monitor,&user,0);CHKERRQ(ierr);
ierr = SNESSetDM(snes,user.packer);CHKERRQ(ierr);
ierr = SNESSolve(snes,NULL,U);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
/* Output the data */
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,"wu_testtwogroups",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
ierr = PetscViewerPushFormat(viewer,PETSC_VIEWER_BINARY_MATLAB);CHKERRQ(ierr);
ierr = DMCompositeGetLocalVectors(user.packer,&vlambda,&vphi1,&vphi2);CHKERRQ(ierr);
ierr = DMCompositeScatter(user.packer,U,vlambda,vphi1,vphi2);CHKERRQ(ierr);
ierr = VecView(vlambda,viewer);CHKERRQ(ierr);
ierr = VecView(vphi1,viewer);CHKERRQ(ierr);
ierr = VecView(vphi2,viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = DMDestroy(&user.red1);CHKERRQ(ierr);
ierr = DMDestroy(&user.da1);CHKERRQ(ierr);
ierr = DMDestroy(&user.da2);CHKERRQ(ierr);
ierr = DMDestroy(&user.packer);CHKERRQ(ierr);
ierr = VecDestroy(&U);CHKERRQ(ierr);
ierr = VecDestroy(&FU);CHKERRQ(ierr);
ierr = VecDestroy(&vlambda);CHKERRQ(ierr);
ierr = VecDestroy(&vphi1);CHKERRQ(ierr);
ierr = VecDestroy(&vphi2);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
PetscErrorCode FormInitialGuess(UserCtx *user,Vec U,Mat A1,Mat A2)
{
PetscErrorCode ierr;
PetscInt xs,xm,i,N;
PetscInt ys,ym,j,M;/*y direction index*/
PetscInt id;
PetscScalar *lambda,**phi1,**phi2;
PetscScalar *initial_phi1,*initial_phi2;
Vec vlambda,vphi1,vphi2;
PetscFunctionBeginUser;
ierr = MatSeqAIJGetArray(A1,&initial_phi1);CHKERRQ(ierr);
ierr = MatSeqAIJGetArray(A2,&initial_phi2);CHKERRQ(ierr);
ierr = DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr = DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = DMDAGetCorners(user->da1,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);//get two direction index
ierr = DMDAGetInfo(user->da1,0,&N,&M,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
if (xs == 0 && ys == 0) lambda[0] = 1.0;
for ( j=ys ; j < ys+ym ; j++ ){
for ( i=xs ; i < xs+xm ; i++ ){
id = j*xm + i ;
phi1[j][i] = initial_phi1[id] ;
phi2[j][i] = initial_phi2[id] ;
}
}
ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = DMCompositeGather(user->packer,INSERT_VALUES,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = DMCompositeRestoreLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
/*
Evaluates FU = Gradiant(L(w,u,lambda))
*/
PetscErrorCode FormFunction(SNES snes,Vec U,Vec FU,void *dummy)
{
UserCtx *user = (UserCtx*)dummy;
PetscErrorCode ierr;
PetscInt xs,xm,i,N;
PetscInt ys,ym,j,M;/*y direction index*/
PetscInt xints,xinte,yints,yinte;
PetscScalar *lambda,**phi1,**phi2,*flambda,**fphi1,**fphi2;
PetscScalar dx,dy;
PetscScalar dotphi;//wu-define
PetscScalar rightside,e,b,c,dd,a;
PetscReal phi1boundary,phi2boundary;
Vec vlambda,vphi1,vphi2,vflambda,vfphi1,vfphi2;
PetscFunctionBeginUser;
ierr = DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr = DMCompositeGetLocalVectors(user->packer,&vfphi1,&vfphi2,&vflambda);CHKERRQ(ierr);
ierr = DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = DMDAGetCorners(user->da1,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);//get two direction index
ierr = DMDAGetInfo(user->da1,0,&N,&M,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = VecGetArray(vflambda,&flambda);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da1,vfphi1,&fphi1);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da2,vfphi2,&fphi2);CHKERRQ(ierr);
dx = 5.0;
dy = 5.0;
phi1boundary = 0.0;
phi2boundary = 0.0;
/* residual f_lambda */
if (xs == 0 && ys == 0 ) { /* only first processor computes this */
dotphi = 0.0;
for ( j=ys ; j < ys+ym ; j++ ){
for ( i=xs ; i < xs+xm ; i++ ){
dotphi += phi1[j][i]*phi1[j][i] + phi2[j][i]*phi2[j][i] ;
}
}
flambda[0] = dotphi - 1.0 ;
}
/* bottom boundary */
if (ys == 0) {
j = 0;
yints = ys + 1;
/* bottom edge */
for (i=xs; i<xs+xm; i++) {
/*phi1 field*/
rightside = phi1boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi1[j][i] = a * phi1[j][i] - rightside;
/*phi2 field*/
rightside = phi2boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi2[j][i] = a * phi2[j][i] - rightside;
}
}
/* top boundary */
if (yints == 1) {
j = ys + ym - 1;
yinte = ys + ym - 1;
for (i=xs; i<xs+xm; i++) {
/*phi1 field*/
rightside = phi1boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi1[j][i] = a * phi1[j][i] - rightside;
/*phi2 field*/
rightside = phi2boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi2[j][i] = a * phi2[j][i] - rightside;
}
}
/* left boundary */
if (xs == 0) {
i = 0;
xints = xs + 1;
/* bottom edge */
for (j=yints; j<yinte; j++) {
/*phi1 field*/
rightside = phi1boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi1[j][i] = a * phi1[j][i] - rightside;
/*phi2 field*/
rightside = phi2boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi2[j][i] = a * phi2[j][i] - rightside;
}
}
/* right boundary */
if (xints == 1) {
i = xs + xm - 1;
xinte = xs + xm - 1;
for (j=yints; j<yinte; j++) {
/*phi1 field*/
rightside = phi1boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi1[j][i] = a * phi1[j][i] - rightside;
/*phi2 field*/
rightside = phi2boundary;
b = 0.0;
c = 0.0;
dd = 0.0;
e = 0.0;
a = 1.0;
fphi2[j][i] = a * phi2[j][i] - rightside;
}
}
/* inner points */
for (j=yints; j<yinte; j++){
for (i=xints; i<xinte; i++){
/*phi1 field*/
rightside = - (0.0121 + 0.0241) * phi1[j][i] + lambda[0] * 0.0085 * phi1[j][i] + lambda[0] * 0.1851 * phi2[j][i] ;
b = 1.0 / ( dx * dx / (2* 1.267) + dx * dx / (2* 1.267 ) ) ;
c = 1.0 / ( dx * dx / (2* 1.267) + dx * dx / (2* 1.267 ) ) ;
dd = 1.0/ ( dy * dy / (2* 1.267) + dy * dy / (2* 1.267 ) ) ;
e = 1.0 / ( dy * dy / (2* 1.267) + dy * dy / (2* 1.267 ) ) ;
a = b + c + dd + e ;
fphi1[j][i] = a* phi1[j][i] - b* phi1[j][i-1] - c* phi1[j][i+1] - dd* phi1[j-1][i] - e* phi1[j+1][i] - rightside ;
/*phi2 field*/
rightside = - 0.121 * phi2[j][i] + 0.0241 * phi1[j][i] ;
b = 1.0 / ( dx * dx / (2* 0.354) + dx * dx / (2* 0.354 ) ) ;
c = 1.0 / ( dx * dx / (2* 0.354) + dx * dx / (2* 0.354 ) ) ;
dd = 1.0/ ( dy * dy / (2* 0.354) + dy * dy / (2* 0.354 ) ) ;
e = 1.0 / ( dy * dy / (2* 0.354) + dy * dy / (2* 0.354 ) ) ;
a = b + c + dd + e ;
fphi2[j][i] = a* phi2[j][i] - b* phi2[j][i-1] - c* phi2[j][i+1] - dd* phi2[j-1][i] - e* phi2[j+1][i] - rightside ;
}
}
ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = VecRestoreArray(vflambda,&flambda);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da1,vfphi1,&fphi1);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da2,vfphi2,&fphi2);CHKERRQ(ierr);
ierr = DMCompositeGather(user->packer,INSERT_VALUES,FU,vfphi1,vfphi2,vflambda);CHKERRQ(ierr);
ierr = DMCompositeRestoreLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr = DMCompositeRestoreLocalVectors(user->packer,&vfphi1,&vfphi2,&vflambda);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
PetscErrorCode FormJacobian(SNES snes,Vec U,Mat J,Mat B,void *ctx)
{
UserCtx *user = (UserCtx*)ctx;
PetscErrorCode ierr;
PetscInt xs,xm,i,N;
PetscInt ys,ym,j,M;
PetscInt row1,col1;
MatStencil col[5],row;
Mat Bll,B11,B22,B12,B21;
IS *is;
PetscScalar *lambda,**phi1,**phi2;
PetscScalar unit;
Vec vlambda,vphi1,vphi2;
PetscFunctionBeginUser;
// unit = 1.0;
// row1 = 0;
// for (int i = 0; i < 961; i++)
// {
// ierr = MatSetValues(B,1,&row1,1,&row1,&unit,INSERT_VALUES);CHKERRQ(ierr);
// row1++;
// }
ierr = DMDAGetCorners(user->da1,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);//get two direction index
ierr = DMDAGetInfo(user->da1,0,&N,&M,0,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr);
ierr = DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr = DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecGetArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = DMCompositeGetLocalISs(user->packer,&is);CHKERRQ(ierr);
ierr = MatGetLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
ierr = MatGetLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
// ierr = MatGetLocalSubMatrix(B,is[2],is[2],&Bll);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.j = j; row.i = i;
unit = 1.0;
ierr = MatSetValuesStencil(B11,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
}
}
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
row.j = j; row.i = i;
unit = 1.0;
ierr = MatSetValuesStencil(B22,1,&row,1,&row,&unit,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatRestoreLocalSubMatrix(B,is[0],is[0],&B11);CHKERRQ(ierr);
ierr = MatRestoreLocalSubMatrix(B,is[1],is[1],&B22);CHKERRQ(ierr);
unit = 1.0;
row1 = 960;//last row global index
ierr = MatSetValues(B,1,&row1,1,&row1,&unit,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da1,vphi1,&phi1);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(user->da2,vphi2,&phi2);CHKERRQ(ierr);
ierr = ISDestroy(&is[0]);CHKERRQ(ierr);
ierr = ISDestroy(&is[1]);CHKERRQ(ierr);
ierr = ISDestroy(&is[2]);CHKERRQ(ierr);
ierr = PetscFree(is);CHKERRQ(ierr);
ierr = DMCompositeRestoreLocalVectors(user->packer,&vlambda,&vphi1,&vphi2);CHKERRQ(ierr);
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd (B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (J != B) {
ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd (J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *dummy)
{
UserCtx *user = (UserCtx*)dummy;
PetscErrorCode ierr;
Vec vlambda,vphi1,vphi2,U;
PetscReal *lambda;
PetscReal keff;
PetscFunctionBeginUser;
ierr = PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);CHKERRQ(ierr);
ierr = SNESGetSolution(snes,&U);CHKERRQ(ierr);
ierr = DMCompositeGetLocalVectors(user->packer,&vphi1,&vphi2,&vlambda);CHKERRQ(ierr);
ierr = DMCompositeScatter(user->packer,U,vphi1,vphi2,vlambda);CHKERRQ(ierr);
ierr = VecGetArray(vlambda,&lambda);CHKERRQ(ierr);
keff = 1.0 / lambda[0] ;
ierr = VecRestoreArray(vlambda,&lambda);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"iter = %D, Keff ======= %g\n",its,(double)keff);CHKERRQ(ierr);
//ierr = VecView(vphi1,PETSC_VIEWER_DRAW_WORLD);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: log
Type: application/octet-stream
Size: 2556 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181105/59aceb11/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: wu-readtwogroups
Type: application/octet-stream
Size: 12000 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20181105/59aceb11/attachment-0003.obj>
More information about the petsc-users
mailing list