本文整理汇总了C++中VecView函数的典型用法代码示例。如果您正苦于以下问题：C++ VecView函数的具体用法？C++ VecView怎么用？C++ VecView使用的例子？那么, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了VecView函数的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: test_axpy_dot_max
PetscErrorCode test_axpy_dot_max( void )
{
Vec x1,y1, x2,y2;
Vec tmp_buf[2];
Vec X, Y;
PetscReal real,real2;
PetscScalar scalar;
PetscInt index;
PetscErrorCode ierr;
PetscFunctionBegin;
PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME );
gen_test_vector( PETSC_COMM_WORLD, 4, 0, 1, &x1 );
gen_test_vector( PETSC_COMM_WORLD, 5, 10, 2, &x2 );
gen_test_vector( PETSC_COMM_WORLD, 4, 4, 3, &y1 );
gen_test_vector( PETSC_COMM_WORLD, 5, 5, 1, &y2 );
tmp_buf[0] = x1;
tmp_buf[1] = x2;
ierr = VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&X);CHKERRQ(ierr);
ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
ierr = VecAssemblyEnd(X);CHKERRQ(ierr);
ierr = VecDestroy(&x1);CHKERRQ(ierr);
ierr = VecDestroy(&x2);CHKERRQ(ierr);
tmp_buf[0] = y1;
tmp_buf[1] = y2;
ierr = VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&Y);CHKERRQ(ierr);
ierr = VecAssemblyBegin(Y);CHKERRQ(ierr);
ierr = VecAssemblyEnd(Y);CHKERRQ(ierr);
ierr = VecDestroy(&y1);CHKERRQ(ierr);
ierr = VecDestroy(&y2);CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "VecAXPY \n");
ierr = VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
ierr = VecNestGetSubVec( Y, 0, &y1 );CHKERRQ(ierr);
ierr = VecNestGetSubVec( Y, 1, &y2 );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(1) y1 = \n" );
ierr = VecView( y1, PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(1) y2 = \n" );
ierr = VecView( y2, PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);
ierr = VecDot( X,Y, &scalar );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
ierr = VecDotNorm2( X,Y, &scalar, &real2 );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), real2);
ierr = VecAXPY( Y, 1.0, X ); /* Y <- a X + Y */
ierr = VecNestGetSubVec( Y, 0, &y1 );CHKERRQ(ierr);
ierr = VecNestGetSubVec( Y, 1, &y2 );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(2) y1 = \n" );
ierr = VecView( y1, PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(2) y2 = \n" );
ierr = VecView( y2, PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);
ierr = VecDot( X,Y, &scalar );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi \n", PetscRealPart(scalar), PetscImaginaryPart(scalar) );
ierr = VecDotNorm2( X,Y, &scalar, &real2 );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "X.Y = %lf + %lfi norm2(Y) = %lf\n", PetscRealPart(scalar), PetscImaginaryPart(scalar), real2);
ierr = VecMax( X, &index, &real );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", real, index );
ierr = VecMin( X, &index, &real );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", real, index );
ierr = VecDestroy(&X);CHKERRQ(ierr);
ierr = VecDestroy(&Y);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例2: main
int main(int argc,char **args)
{
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt N0=50,N1=20,N=N0*N1,DIM;
PetscRandom rdm;
PetscScalar a;
PetscReal enorm;
Vec x,y,z;
PetscBool view=PETSC_FALSE,use_interface=PETSC_TRUE;
ierr = PetscInitialize(&argc,&args,(char*)0,help);CHKERRQ(ierr);
#if !defined(PETSC_USE_COMPLEX)
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires complex numbers");
#endif
ierr = PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "FFTW Options", "ex143");CHKERRQ(ierr);
ierr = PetscOptionsBool("-vec_view draw", "View the vectors", "ex143", view, &view, NULL);CHKERRQ(ierr);
ierr = PetscOptionsBool("-use_FFTW_interface", "Use PETSc-FFTW interface", "ex143",use_interface, &use_interface, NULL);CHKERRQ(ierr);
ierr = PetscOptionsEnd();CHKERRQ(ierr);
ierr = PetscOptionsGetBool(NULL,"-use_FFTW_interface",&use_interface,NULL);CHKERRQ(ierr);
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
if (!use_interface) {
/* Use mpi FFTW without PETSc-FFTW interface, 2D case only */
/*---------------------------------------------------------*/
fftw_plan fplan,bplan;
fftw_complex *data_in,*data_out,*data_out2;
ptrdiff_t alloc_local,local_n0,local_0_start;
DIM = 2;
if (!rank) {
ierr = PetscPrintf(PETSC_COMM_SELF,"Use FFTW without PETSc-FFTW interface, DIM %D\n",DIM);CHKERRQ(ierr);
}
fftw_mpi_init();
N = N0*N1;
alloc_local = fftw_mpi_local_size_2d(N0,N1,PETSC_COMM_WORLD,&local_n0,&local_0_start);
data_in = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
data_out = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
data_out2 = (fftw_complex*)fftw_malloc(sizeof(fftw_complex)*alloc_local);
ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_in,&x);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) x, "Real Space vector");CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out,&y);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr);
ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,(PetscInt)local_n0*N1,(PetscInt)N,(const PetscScalar*)data_out2,&z);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr);
fplan = fftw_mpi_plan_dft_2d(N0,N1,data_in,data_out,PETSC_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE);
bplan = fftw_mpi_plan_dft_2d(N0,N1,data_out,data_out2,PETSC_COMM_WORLD,FFTW_BACKWARD,FFTW_ESTIMATE);
ierr = VecSetRandom(x, rdm);CHKERRQ(ierr);
if (view) {ierr = VecView(x,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
fftw_execute(fplan);
if (view) {ierr = VecView(y,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
fftw_execute(bplan);
/* Compare x and z. FFTW computes an unnormalized DFT, thus z = N*x */
a = 1.0/(PetscReal)N;
ierr = VecScale(z,a);CHKERRQ(ierr);
if (view) {ierr = VecView(z, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
ierr = VecAXPY(z,-1.0,x);CHKERRQ(ierr);
ierr = VecNorm(z,NORM_1,&enorm);CHKERRQ(ierr);
if (enorm > 1.e-11 && !rank) {
ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %g\n",(double)enorm);CHKERRQ(ierr);
}
/* Free spaces */
fftw_destroy_plan(fplan);
fftw_destroy_plan(bplan);
fftw_free(data_in); ierr = VecDestroy(&x);CHKERRQ(ierr);
fftw_free(data_out); ierr = VecDestroy(&y);CHKERRQ(ierr);
fftw_free(data_out2);ierr = VecDestroy(&z);CHKERRQ(ierr);
} else {
/* Use PETSc-FFTW interface */
/*-------------------------------------------*/
PetscInt i,*dim,k;
Mat A;
N=1;
for (i=1; i<5; i++) {
DIM = i;
ierr = PetscMalloc1(i,&dim);CHKERRQ(ierr);
for (k=0; k<i; k++) {
dim[k]=30;
}
N *= dim[i-1];
/* Create FFTW object */
if (!rank) printf("Use PETSc-FFTW interface...%d-DIM: %d\n",(int)DIM,(int)N);
//.........这里部分代码省略.........
示例3: main
//.........这里部分代码省略.........
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = PetscLogStagePop();CHKERRQ(ierr);
/*
Create parallel vectors compatible with the DMDA.
*/
ierr = DMCreateGlobalVector(da,&u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
/*
Set exact solution; then compute right-hand-side vector.
By default we use an exact solution of a vector with all
elements of 1.0; Alternatively, using the runtime option
-random_sol forms a solution vector with random components.
*/
ierr = PetscOptionsGetBool(NULL,"-random_exact_sol",&flg,NULL);CHKERRQ(ierr);
if (flg) {
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
} else {
ierr = VecSet(u,1.);CHKERRQ(ierr);
}
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/*
View the exact solution vector if desired
*/
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(NULL,"-view_exact_sol",&flg,NULL);CHKERRQ(ierr);
if (flg) {ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -示例4: efs_solve
PetscErrorCode efs_solve(efs *slv)
{
PetscErrorCode ierr;
PetscInt i, j, k, xs, ys, zs, xm, ym, zm;
PetscInt ngx,ngy,ngz;
DMBoundaryType bx, by, bz;
slv->ts++;
ierr = KSPSolve(slv->ksp, NULL, NULL);CHKERRQ(ierr);
if (efs_log(slv, EFS_LOG_STATUS)) {
ierr = ef_io_print(slv->comm, "Solve completed");CHKERRQ(ierr);
}
ierr = KSPGetSolution(slv->ksp, &slv->x);CHKERRQ(ierr);
ierr = DMDAGetInfo(slv->dm, 0, &ngx, &ngy, &ngz,0,0,0,0,0, &bx, &by, &bz,0);CHKERRQ(ierr);
if (efs_log(slv, EFS_LOG_PROBLEM)) {
PetscViewer xout;
ierr = PetscViewerASCIIOpen(slv->comm, "x.txt", &xout);CHKERRQ(ierr);
ierr = VecView(slv->x, xout);CHKERRQ(ierr);
}
ierr = DMDAGetCorners(slv->dm, &xs, &ys, &zs, &xm, &ym, &zm);CHKERRQ(ierr);
if (slv->grid.nd == 2) {
Vec loc_x;
PetscScalar **xvec;
double **phi;
ierr = ef_dmap_get(slv->dmap, slv->state.phi, &phi);CHKERRQ(ierr);
ierr = DMGetLocalVector(slv->dm, &loc_x);CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(slv->dm, slv->x, INSERT_VALUES, loc_x);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(slv->dm, slv->x, INSERT_VALUES, loc_x);CHKERRQ(ierr);
ierr = DMDAVecGetArray(slv->dm, loc_x, &xvec);CHKERRQ(ierr);
if ((by == DM_BOUNDARY_PERIODIC) && (ys + ym == ngy)) {
for (i = xs; i < xs + xm; i++) {
xvec[ngy-1][i] = xvec[ngy][i];
}
}
if ((bx == DM_BOUNDARY_PERIODIC) && (xs + xm == ngx)) {
for (j = ys; j < ys + ym; j++) {
xvec[j][ngx-1] = xvec[j][ngx];
}
}
for (j = ys; j < ys + ym; j++) {
for (i = xs; i < xs + xm; i++) {
phi[j][i] = xvec[j][i];
}
}
ierr = ef_dmap_restore(slv->dmap, &phi);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(slv->dm, loc_x, &xvec);CHKERRQ(ierr);
if (efs_log(slv, EFS_LOG_VTK)) {
ierr = DMLocalToGlobalBegin(slv->dm, loc_x, INSERT_VALUES, slv->x);CHKERRQ(ierr);
ierr = DMLocalToGlobalEnd(slv->dm, loc_x, INSERT_VALUES, slv->x);CHKERRQ(ierr);
ierr = ef_io_vtkwrite(slv->dm, slv->x, "phi", slv->grid.id, slv->ts);CHKERRQ(ierr);
}
} else if (slv->grid.nd == 3) {
PetscScalar ***xvec;
double ***phi;
Vec loc_x;
ierr = ef_dmap_get(slv->dmap, slv->state.phi, &phi);CHKERRQ(ierr);
ierr = DMGetLocalVector(slv->dm, &loc_x);CHKERRQ(ierr);
ierr = DMGlobalToLocalBegin(slv->dm, slv->x, INSERT_VALUES, loc_x);CHKERRQ(ierr);
ierr = DMGlobalToLocalEnd(slv->dm, slv->x, INSERT_VALUES, loc_x);CHKERRQ(ierr);
ierr = DMDAVecGetArray(slv->dm, loc_x, &xvec);CHKERRQ(ierr);
if ((bz == DM_BOUNDARY_PERIODIC) && (zs + zm == ngz)) {
for (j = ys; j < ys + ym; j++) {
for (i = xs; i < xs + xm; i++) {
xvec[ngz-1][j][i] = xvec[ngz][j][i];
}
}
}
if ((by == DM_BOUNDARY_PERIODIC) && (ys + ym == ngy)) {
for (k = zs; k < zs + zm; k++) {
for (i = xs; i < xs + xm; i++) {
xvec[k][ngy-1][i] = xvec[k][ngy][i];
}
}
}
if ((bx == DM_BOUNDARY_PERIODIC) && (xs + xm == ngx)) {
for (k = zs; k < zs + zm; k++) {
for (j = ys; j < ys + ym; j++) {
xvec[k][j][ngx-1] = xvec[k][j][ngx];
}
}
}
for (k = zs; k < zs + zm; k++) {
for (j = ys; j < ys + ym; j++) {
for (i = xs; i < xs + xm; i++) {
phi[k][j][i] = xvec[k][j][i];
}
}
}
ierr = ef_dmap_restore(slv->dmap, &phi);CHKERRQ(ierr);
ierr = DMDAVecRestoreArray(slv->dm, loc_x, &xvec);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例5: main
int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* KSP context */
PetscErrorCode ierr;
PetscInt n = 10,its, dim,p = 1,use_random;
PetscScalar none = -1.0,pfive = 0.5;
PetscReal norm;
PetscRandom rctx;
TestType type;
PetscBool flg;
PetscInitialize(&argc,&args,(char *)0,help);
ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(PETSC_NULL,"-p",&p,PETSC_NULL);CHKERRQ(ierr);
switch (p) {
case 1: type = TEST_1; dim = n; break;
case 2: type = TEST_2; dim = n; break;
case 3: type = TEST_3; dim = n; break;
case 4: type = HELMHOLTZ_1; dim = n*n; break;
case 5: type = HELMHOLTZ_2; dim = n*n; break;
default: type = TEST_1; dim = n;
}
/* Create vectors */
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);
ierr = VecSetSizes(x,PETSC_DECIDE,dim);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&b);CHKERRQ(ierr);
ierr = VecDuplicate(x,&u);CHKERRQ(ierr);
use_random = 1;
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-norandom",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
use_random = 0;
ierr = VecSet(u,pfive);CHKERRQ(ierr);
} else {
ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rctx);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
ierr = VecSetRandom(u,rctx);CHKERRQ(ierr);
}
/* Create and assemble matrix */
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
ierr = FormTestMatrix(A,n,type);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
flg = PETSC_FALSE;
ierr = PetscOptionsGetBool(PETSC_NULL,"-printout",&flg,PETSC_NULL);CHKERRQ(ierr);
if (flg) {
ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(b,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/* Create KSP context; set operators and options; solve linear system */
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
/* Check error */
ierr = VecAXPY(x,none,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %G,Iterations %D\n",norm,its);CHKERRQ(ierr);
/* Free work space */
ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
if (use_random) {ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);}
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例6: main
int main(int argc, char **argv)
{
MPI_Comm comm;
PetscMPIInt rank;
PetscErrorCode ierr;
User user;
PetscLogDouble v1, v2;
PetscInt nplot = 0;
char filename1[2048], fileName[2048];
PetscBool set = PETSC_FALSE;
PetscInt steps_output;
ierr = PetscInitialize(&argc, &argv, (char*) 0, help);CHKERRQ(ierr);
comm = PETSC_COMM_WORLD;
ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
ierr = PetscNew(&user);CHKERRQ(ierr);
ierr = PetscNew(&user->algebra);CHKERRQ(ierr);
ierr = PetscNew(&user->model);CHKERRQ(ierr);
ierr = PetscNew(&user->model->physics);CHKERRQ(ierr);
Algebra algebra = user->algebra;
ierr = LoadOptions(comm, user);CHKERRQ(ierr);
ierr = PetscTime(&v1);CHKERRQ(ierr);
ierr = CreateMesh(comm, user);CHKERRQ(ierr);
ierr = PetscTime(&v2);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,
"Read and Distribute mesh takes %f sec \n", v2 - v1);CHKERRQ(ierr);
ierr = SetUpLocalSpace(user);CHKERRQ(ierr); //Set up the dofs of each element
ierr = ConstructGeometryFVM(&user->facegeom, &user->cellgeom, user);CHKERRQ(ierr);
ierr = LimiterSetup(user);CHKERRQ(ierr);
if(user->output_solution){
// the output file options
ierr = PetscOptionsBegin(PETSC_COMM_WORLD,0,"Options for output solution",0);CHKERRQ(ierr);
ierr = PetscOptionsString("-solutionfile", "solution file", "AeroSim.c", filename1,filename1, 2048, &set);CHKERRQ(ierr);
if(!set){SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,"please use option -solutionfile to specify solution file name \n");}
ierr = PetscOptionsInt("-steps_output", "the number of time steps between two outputs", "", steps_output, &steps_output, &set);CHKERRQ(ierr);
if(!set){ steps_output = 1;}
ierr = PetscOptionsEnd();CHKERRQ(ierr);
}
if (user->TimeIntegralMethod == EXPLICITMETHOD) {
if(user->myownexplicitmethod){
ierr = PetscPrintf(PETSC_COMM_WORLD,"Using the fully explicit method based on my own routing\n");CHKERRQ(ierr);
user->current_time = user->initial_time;
user->current_step = 1;
ierr = DMCreateGlobalVector(user->dm, &algebra->solution);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) algebra->solution, "solution");CHKERRQ(ierr);
ierr = SetInitialCondition(user->dm, algebra->solution, user);CHKERRQ(ierr);
ierr = VecDuplicate(algebra->solution, &algebra->fn);CHKERRQ(ierr);
ierr = VecDuplicate(algebra->solution, &algebra->oldsolution);CHKERRQ(ierr);
if(user->Explicit_RK2){
ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the second order Runge Kutta method \n");CHKERRQ(ierr);
}else{
ierr = PetscPrintf(PETSC_COMM_WORLD,"Use the first order forward Euler method \n");CHKERRQ(ierr);
}
nplot = 0; //the plot step
while(user->current_time < (user->final_time - 0.05 * user->dt)){
user->current_time = user->current_time + user->dt;
ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);
PetscReal fnnorm;
ierr = VecNorm(algebra->fn,NORM_INFINITY,&fnnorm);CHKERRQ(ierr);
if(0){
PetscViewer viewer;
ierr = OutputVTK(user->dm, "function.vtk", &viewer);CHKERRQ(ierr);
ierr = VecView(algebra->fn, viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with founction norm = %g \n",
user->current_step, user->current_time, fnnorm);CHKERRQ(ierr);
//break;
}
if(user->Explicit_RK2){
ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);//U^n
ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);//U^{(1)}
ierr = FormTimeStepFunction(user, algebra, algebra->solution, algebra->fn);CHKERRQ(ierr);//f(U^{(1)})
ierr = VecAXPY(algebra->solution, 1.0, algebra->oldsolution);CHKERRQ(ierr);//U^n + U^{(1)}
ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);// + dt*f(U^{(1)})
ierr = VecScale(algebra->solution, 0.5);CHKERRQ(ierr);
}else{
ierr = VecCopy(algebra->solution, algebra->oldsolution);CHKERRQ(ierr);
ierr = VecAXPY(algebra->solution, user->dt, algebra->fn);CHKERRQ(ierr);
}
{// Monitor the solution and function norms
PetscReal norm;
PetscLogDouble space =0;
PetscInt size;
ierr = VecNorm(algebra->solution,NORM_INFINITY,&norm);CHKERRQ(ierr);
ierr = VecGetSize(algebra->solution, &size);CHKERRQ(ierr);
norm = norm/size;
if (norm>1.e5) {
SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_LIB,
"The norm of the solution is: %f (current time: %f). The explicit method is going to DIVERGE!!!", norm, user->current_time);
}
if (user->current_step%10==0) {
ierr = PetscPrintf(PETSC_COMM_WORLD,"Step %D at time %g with solution norm = %g and founction norm = %g \n",
//.........这里部分代码省略.........
示例7: main
PetscInt main(PetscInt argc,char **args)
{
PetscErrorCode ierr;
PetscMPIInt rank,size;
PetscInt N0=3,N1=3,N2=3,N=N0*N1*N2;
PetscRandom rdm;
PetscScalar a;
PetscReal enorm;
Vec x,y,z,input,output;
PetscBool view=PETSC_FALSE,use_interface=PETSC_TRUE;
Mat A;
PetscInt DIM, dim[3],vsize;
PetscReal fac;
ierr = PetscInitialize(&argc,&args,(char *)0,help);CHKERRQ(ierr);
#if defined(PETSC_USE_COMPLEX)
SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP, "This example requires real numbers");
#endif
ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size);CHKERRQ(ierr);
ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
ierr = PetscRandomCreate(PETSC_COMM_WORLD, &rdm);CHKERRQ(ierr);
ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&input);CHKERRQ(ierr);
ierr = VecSetSizes(input,PETSC_DECIDE,N);CHKERRQ(ierr);
ierr = VecSetFromOptions(input);CHKERRQ(ierr);
ierr = VecSetRandom(input,rdm);CHKERRQ(ierr);
ierr = VecDuplicate(input,&output);
// ierr = VecGetSize(input,&vsize);CHKERRQ(ierr);
// printf("Size of the input Vector is %d\n",vsize);
DIM = 3;
dim[0] = N0; dim[1] = N1; dim[2] = N2;
ierr = MatCreateFFT(PETSC_COMM_WORLD,DIM,dim,MATFFTW,&A);CHKERRQ(ierr);
ierr = MatGetVecs(A,&x,&y);CHKERRQ(ierr);
ierr = MatGetVecs(A,&z,PETSC_NULL);CHKERRQ(ierr);
ierr = VecGetSize(y,&vsize);CHKERRQ(ierr);
printf("The vector size from the main routine is %d\n",vsize);
ierr = InputTransformFFT(A,input,x);CHKERRQ(ierr);
ierr = MatMult(A,x,y);CHKERRQ(ierr);
ierr = MatMultTranspose(A,y,z);CHKERRQ(ierr);
ierr = OutputTransformFFT(A,z,output);CHKERRQ(ierr);
fac = 1.0/(PetscReal)N;
ierr = VecScale(output,fac);CHKERRQ(ierr);
ierr = VecAssemblyBegin(input);CHKERRQ(ierr);
ierr = VecAssemblyEnd(input);CHKERRQ(ierr);
ierr = VecAssemblyBegin(output);CHKERRQ(ierr);
ierr = VecAssemblyEnd(output);CHKERRQ(ierr);
ierr = VecView(input,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(output,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecAXPY(output,-1.0,input);CHKERRQ(ierr);
ierr = VecNorm(output,NORM_1,&enorm);CHKERRQ(ierr);
// if (enorm > 1.e-14){
if (!rank)
ierr = PetscPrintf(PETSC_COMM_SELF," Error norm of |x - z| %e\n",enorm);CHKERRQ(ierr);
// }
// ierr = MatGetVecs(A,&z,PETSC_NULL);CHKERRQ(ierr);
// printf("Vector size from ex148 %d\n",vsize);
// ierr = PetscObjectSetName((PetscObject) x, "Real space vector");CHKERRQ(ierr);
// ierr = PetscObjectSetName((PetscObject) y, "Frequency space vector");CHKERRQ(ierr);
// ierr = PetscObjectSetName((PetscObject) z, "Reconstructed vector");CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}示例8: main
//.........这里部分代码省略.........
ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
ierr = TSSetType(ts,TSRK);CHKERRQ(ierr);
ierr = TSSetRHSFunction(ts,NULL,(TSRHSFunction)RHSFunction,&ctx);CHKERRQ(ierr);
ierr = TSSetRHSJacobian(ts,A,A,(TSRHSJacobian)RHSJacobian,&ctx);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set initial conditions
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetSolution(ts,U);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Save trajectory of solution so that TSAdjointSolve() may be used
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
ierr = MatCreateVecs(A,&lambda[0],NULL);CHKERRQ(ierr);
/* Set initial conditions for the adjoint integration */
ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr);
y_ptr[0] = 0.0; y_ptr[1] = 0.0;
ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr);
ierr = MatCreateVecs(Jacp,&mu[0],NULL);CHKERRQ(ierr);
ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = -1.0;
ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr);
ierr = TSSetCostGradients(ts,1,lambda,mu);CHKERRQ(ierr);
ierr = TSSetCostIntegrand(ts,1,(PetscErrorCode (*)(TS,PetscReal,Vec,Vec,void*))CostIntegrand,
(PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDYFunction,
(PetscErrorCode (*)(TS,PetscReal,Vec,Vec*,void*))DRDPFunction,PETSC_TRUE,&ctx);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set solver options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = TSSetDuration(ts,PETSC_DEFAULT,10.0);CHKERRQ(ierr);
ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,0.0,.01);CHKERRQ(ierr);
ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve nonlinear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
if (ensemble) {
for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
ierr = VecGetArray(U,&u);CHKERRQ(ierr);
u[0] = PetscAsinScalar(ctx.Pm/ctx.Pmax);
u[1] = ctx.omega_s;
u[0] += du[0];
u[1] += du[1];
ierr = VecRestoreArray(U,&u);CHKERRQ(ierr);
ierr = TSSetInitialTimeStep(ts,0.0,.01);CHKERRQ(ierr);
ierr = TSSolve(ts,U);CHKERRQ(ierr);
}
} else {
ierr = TSSolve(ts,U);CHKERRQ(ierr);
}
ierr = VecView(U,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Adjoint model starts here
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/* Set initial conditions for the adjoint integration */
ierr = VecGetArray(lambda[0],&y_ptr);CHKERRQ(ierr);
y_ptr[0] = 0.0; y_ptr[1] = 0.0;
ierr = VecRestoreArray(lambda[0],&y_ptr);CHKERRQ(ierr);
ierr = VecGetArray(mu[0],&x_ptr);CHKERRQ(ierr);
x_ptr[0] = -1.0;
ierr = VecRestoreArray(mu[0],&x_ptr);CHKERRQ(ierr);
/* Set RHS JacobianP */
ierr = TSAdjointSetRHSJacobian(ts,Jacp,RHSJacobianP,&ctx);CHKERRQ(ierr);
ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n");CHKERRQ(ierr);
ierr = VecView(lambda[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecView(mu[0],PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = TSGetCostIntegral(ts,&q);CHKERRQ(ierr);
ierr = VecView(q,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = VecGetArray(q,&x_ptr);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n cost function=%g\n",(double)(x_ptr[0]-ctx.Pm));CHKERRQ(ierr);
ierr = VecRestoreArray(q,&x_ptr);CHKERRQ(ierr);
ierr = ComputeSensiP(lambda[0],mu[0],&ctx);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Free work space. All PETSc objects should be destroyed when they are no longer needed.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = MatDestroy(&A);CHKERRQ(ierr);
ierr = MatDestroy(&Jacp);CHKERRQ(ierr);
ierr = VecDestroy(&U);CHKERRQ(ierr);
ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr);
ierr = VecDestroy(&mu[0]);CHKERRQ(ierr);
ierr = TSDestroy(&ts);CHKERRQ(ierr);
ierr = PetscFinalize();
return(0);
}示例9: main
int main(int argc,char **argv)
{
KSP ksp;
PC pc;
Mat A,M;
Vec X,B,D;
MPI_Comm comm;
PetscScalar v;
KSPConvergedReason reason;
PetscInt i,j,its;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscInitialize(&argc,&argv,0,help);CHKERRQ(ierr);
ierr = PetscOptionsSetValue("-options_left",PETSC_NULL);CHKERRQ(ierr);
comm = MPI_COMM_SELF;
/*
* Construct the Kershaw matrix
* and a suitable rhs / initial guess
*/
ierr = MatCreateSeqAIJ(comm,4,4,4,0,&A);CHKERRQ(ierr);
ierr = VecCreateSeq(comm,4,&B);CHKERRQ(ierr);
ierr = VecDuplicate(B,&X);CHKERRQ(ierr);
for (i=0; i<4; i++) {
v=3;
ierr = MatSetValues(A,1,&i,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
v=1;
ierr = VecSetValues(B,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
}
i=0; v=0;
ierr = VecSetValues(X,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
for (i=0; i<3; i++) {
v=-2; j=i+1;
ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
}
i=0; j=3; v=2;
ierr = MatSetValues(A,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatSetValues(A,1,&j,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = VecAssemblyBegin(B);CHKERRQ(ierr);
ierr = VecAssemblyEnd(B);CHKERRQ(ierr);
printf("\nThe Kershaw matrix:\n\n"); MatView(A,0);
/*
* A Conjugate Gradient method
* with ILU(0) preconditioning
*/
ierr = KSPCreate(comm,&ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
ierr = KSPSetType(ksp,KSPCG);CHKERRQ(ierr);
ierr = KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);CHKERRQ(ierr);
/*
* ILU preconditioner;
* The iterative method will break down unless you comment in the SetShift
* line below, or use the -pc_factor_shift_positive_definite option.
* Run the code twice: once as given to see the negative pivot and the
* divergence behaviour, then comment in the Shift line, or add the
* command line option, and see that the pivots are all positive and
* the method converges.
*/
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCICC);CHKERRQ(ierr);
/* ierr = PCFactorSetShiftType(prec,MAT_SHIFT_POSITIVE_DEFINITE);CHKERRQ(ierr); */
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
ierr = KSPSetUp(ksp);CHKERRQ(ierr);
/*
* Now that the factorisation is done, show the pivots;
* note that the last one is negative. This in itself is not an error,
* but it will make the iterative method diverge.
*/
ierr = PCFactorGetMatrix(pc,&M);CHKERRQ(ierr);
ierr = VecDuplicate(B,&D);CHKERRQ(ierr);
ierr = MatGetDiagonal(M,D);CHKERRQ(ierr);
printf("\nPivots:\n\n"); VecView(D,0);
/*
* Solve the system;
* without the shift this will diverge with
* an indefinite preconditioner
*/
ierr = KSPSolve(ksp,B,X);CHKERRQ(ierr);
ierr = KSPGetConvergedReason(ksp,&reason);CHKERRQ(ierr);
if (reason==KSP_DIVERGED_INDEFINITE_PC) {
printf("\nDivergence because of indefinite preconditioner;\n");
printf("Run the executable again but with -pc_factor_shift_positive_definite option.\n");
} else if (reason<0) {
printf("\nOther kind of divergence: this should not happen.\n");
} else {
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
printf("\nConvergence in %d iterations.\n",(int)its);
//.........这里部分代码省略.........
示例10: private_VecView_Swarm_XDMF
PetscErrorCode private_VecView_Swarm_XDMF(Vec x,PetscViewer viewer)
{
long int *bytes = NULL;
PetscContainer container = NULL;
const char *viewername;
char datafile[PETSC_MAX_PATH_LEN];
PetscViewer fviewer;
PetscInt N,bs;
const char *vecname;
char fieldname[PETSC_MAX_PATH_LEN];
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
if (container) {
ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
} else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr);
ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr);
/* re-open a sub-viewer for all data fields */
/* name is viewer.name + "_swarm_fields.pbin" */
ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr);
ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_APPEND);CHKERRQ(ierr);
ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr);
ierr = VecGetSize(x,&N);CHKERRQ(ierr);
ierr = VecGetBlockSize(x,&bs);CHKERRQ(ierr);
N = N/bs;
ierr = PetscObjectGetName((PetscObject)x,&vecname);CHKERRQ(ierr);
if (!vecname) {
ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"swarmfield_%D",((PetscObject)x)->tag);CHKERRQ(ierr);
} else {
ierr = PetscSNPrintf(fieldname,PETSC_MAX_PATH_LEN-1,"%s",vecname);CHKERRQ(ierr);
}
/* write data header */
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"<Attribute Center=\"Node\" Name=\"%s\" Type=\"None\">\n",fieldname);CHKERRQ(ierr);
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
if (bs == 1) {
ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D\" Seek=\"%D\">\n",N,bytes[0]);CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",N,bs,bytes[0]);CHKERRQ(ierr);
}
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</Attribute>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
/* write data */
ierr = VecView(x,fviewer);CHKERRQ(ierr);
bytes[0] += sizeof(PetscReal) * N * bs;
ierr = PetscViewerDestroy(&fviewer);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例11: private_DMSwarmView_XDMF
PetscErrorCode private_DMSwarmView_XDMF(DM dm,PetscViewer viewer)
{
PetscBool isswarm = PETSC_FALSE;
const char *viewername;
char datafile[PETSC_MAX_PATH_LEN];
PetscViewer fviewer;
PetscInt k,ng,dim;
Vec dvec;
long int *bytes = NULL;
PetscContainer container = NULL;
const char *dmname;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscObjectQuery((PetscObject)viewer,"XDMFViewerContext",(PetscObject*)&container);CHKERRQ(ierr);
if (container) {
ierr = PetscContainerGetPointer(container,(void**)&bytes);CHKERRQ(ierr);
} else SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Valid to find attached data XDMFViewerContext");
ierr = PetscObjectTypeCompare((PetscObject)dm,DMSWARM,&isswarm);CHKERRQ(ierr);
if (!isswarm) SETERRQ(PetscObjectComm((PetscObject)viewer),PETSC_ERR_SUP,"Only valid for DMSwarm");
ierr = PetscObjectCompose((PetscObject)viewer,"DMSwarm",(PetscObject)dm);CHKERRQ(ierr);
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscObjectGetName((PetscObject)dm,&dmname);CHKERRQ(ierr);
if (!dmname) {
ierr = DMGetOptionsPrefix(dm,&dmname);CHKERRQ(ierr);
}
if (!dmname) {
ierr = PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm\" GridType=\"Uniform\">\n");CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer,"<Grid Name=\"DMSwarm[%s]\" GridType=\"Uniform\">\n",dmname);CHKERRQ(ierr);
}
/* create a sub-viewer for topology, geometry and all data fields */
/* name is viewer.name + "_swarm_fields.pbin" */
ierr = PetscViewerCreate(PetscObjectComm((PetscObject)viewer),&fviewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(fviewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
ierr = PetscViewerBinarySetSkipHeader(fviewer,PETSC_TRUE);CHKERRQ(ierr);
ierr = PetscViewerBinarySetSkipInfo(fviewer,PETSC_TRUE);CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(fviewer,FILE_MODE_WRITE);CHKERRQ(ierr);
ierr = PetscViewerFileGetName(viewer,&viewername);CHKERRQ(ierr);
ierr = private_CreateDataFileNameXDMF(viewername,datafile);CHKERRQ(ierr);
ierr = PetscViewerFileSetName(fviewer,datafile);CHKERRQ(ierr);
ierr = DMSwarmGetSize(dm,&ng);CHKERRQ(ierr);
/* write topology header */
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"<Topology Dimensions=\"%D\" TopologyType=\"Mixed\">\n",ng);CHKERRQ(ierr);
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Int\" Dimensions=\"%D\" Seek=\"%D\">\n",ng*3,bytes[0]);CHKERRQ(ierr);
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</Topology>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
/* write topology data */
for (k=0; k<ng; k++) {
PetscInt pvertex[3];
pvertex[0] = 1;
pvertex[1] = 1;
pvertex[2] = k;
ierr = PetscViewerBinaryWrite(fviewer,pvertex,3,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr);
}
bytes[0] += sizeof(PetscInt) * ng * 3;
/* write geometry header */
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
switch (dim) {
case 1:
SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for 1D");
break;
case 2:
ierr = PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XY\">\n");CHKERRQ(ierr);
break;
case 3:
ierr = PetscViewerASCIIPrintf(viewer,"<Geometry Type=\"XYZ\">\n");CHKERRQ(ierr);
break;
}
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"<DataItem Format=\"Binary\" Endian=\"Big\" DataType=\"Float\" Precision=\"8\" Dimensions=\"%D %D\" Seek=\"%D\">\n",ng,dim,bytes[0]);CHKERRQ(ierr);
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"%s\n",datafile);CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</DataItem>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer,"</Geometry>\n");CHKERRQ(ierr);
ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
/* write geometry data */
ierr = DMSwarmCreateGlobalVectorFromField(dm,DMSwarmPICField_coor,&dvec);CHKERRQ(ierr);
ierr = VecView(dvec,fviewer);CHKERRQ(ierr);
//.........这里部分代码省略.........
示例12: main
int main(int argc,char **args)
{
Mat C;
PetscInt i,j,m = 3,n = 3,Ii,J;
PetscErrorCode ierr;
PetscBool flg;
PetscScalar v;
IS perm,iperm;
Vec x,u,b,y;
PetscReal norm,tol=PETSC_SMALL;
MatFactorInfo info;
PetscMPIInt size;
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
if (size != 1) SETERRQ(PETSC_COMM_WORLD,1,"This is a uniprocessor example only!");
ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr);
ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
ierr = MatSetFromOptions(C);CHKERRQ(ierr);
ierr = MatSetUp(C);CHKERRQ(ierr);
ierr = PetscOptionsHasName(NULL,NULL,"-symmetric",&flg);CHKERRQ(ierr);
if (flg) { /* Treat matrix as symmetric only if we set this flag */
ierr = MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr);
ierr = MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr);
}
/* Create the matrix for the five point stencil, YET AGAIN */
for (i=0; i<m; i++) {
for (j=0; j<n; j++) {
v = -1.0; Ii = j + n*i;
if (i>0) {J = Ii - n; ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
v = 4.0; ierr = MatSetValues(C,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatGetOrdering(C,MATORDERINGRCM,&perm,&iperm);CHKERRQ(ierr);
ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = ISView(perm,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
ierr = VecCreateSeq(PETSC_COMM_SELF,m*n,&u);CHKERRQ(ierr);
ierr = VecSet(u,1.0);CHKERRQ(ierr);
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(u,&y);CHKERRQ(ierr);
ierr = MatMult(C,u,b);CHKERRQ(ierr);
ierr = VecCopy(b,y);CHKERRQ(ierr);
ierr = VecScale(y,2.0);CHKERRQ(ierr);
ierr = MatNorm(C,NORM_FROBENIUS,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF,"Frobenius norm of matrix %g\n",(double)norm);CHKERRQ(ierr);
ierr = MatNorm(C,NORM_1,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF,"One norm of matrix %g\n",(double)norm);CHKERRQ(ierr);
ierr = MatNorm(C,NORM_INFINITY,&norm);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_SELF,"Infinity norm of matrix %g\n",(double)norm);CHKERRQ(ierr);
ierr = MatFactorInfoInitialize(&info);CHKERRQ(ierr);
info.fill = 2.0;
info.dtcol = 0.0;
info.zeropivot = 1.e-14;
info.pivotinblocks = 1.0;
ierr = MatLUFactor(C,perm,iperm,&info);CHKERRQ(ierr);
/* Test MatSolve */
ierr = MatSolve(C,b,x);CHKERRQ(ierr);
ierr = VecView(b,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
ierr = VecView(x,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
if (norm > tol) {
ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolve: Norm of error %g\n",(double)norm);CHKERRQ(ierr);
}
/* Test MatSolveAdd */
ierr = MatSolveAdd(C,b,y,x);CHKERRQ(ierr);
ierr = VecAXPY(x,-1.0,y);CHKERRQ(ierr);
ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
if (norm > tol) {
ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolveAdd(): Norm of error %g\n",(double)norm);CHKERRQ(ierr);
}
ierr = ISDestroy(&perm);CHKERRQ(ierr);
ierr = ISDestroy(&iperm);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr);
ierr = VecDestroy(&y);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = MatDestroy(&C);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}示例13: step
/*
Monitor - User-provided routine to monitor the solution computed at
each timestep. This example plots the solution and computes the
error in two different norms.
Input Parameters:
ts - the timestep context
step - the count of the current step (with 0 meaning the
initial condition)
time - the current time
u - the solution at this timestep
ctx - the user-provided context for this monitoring routine.
In this case we use the application context which contains
information about the problem size, workspace and the exact
solution.
*/
PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
{
AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
PetscErrorCode ierr;
PetscReal en2,en2s,enmax;
PetscDraw draw;
/*
We use the default X windows viewer
PETSC_VIEWER_DRAW_(appctx->comm)
that is associated with the current communicator. This saves
the effort of calling PetscViewerDrawOpen() to create the window.
Note that if we wished to plot several items in separate windows we
would create each viewer with PetscViewerDrawOpen() and store them in
the application context, appctx.
PetscReal buffering makes graphics look better.
*/
ierr = PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw);CHKERRQ(ierr);
ierr = PetscDrawSetDoubleBuffer(draw);CHKERRQ(ierr);
ierr = VecView(u,PETSC_VIEWER_DRAW_(appctx->comm));CHKERRQ(ierr);
/*
Compute the exact solution at this timestep
*/
ierr = ExactSolution(time,appctx->solution,appctx);CHKERRQ(ierr);
/*
Print debugging information if desired
*/
if (appctx->debug) {
ierr = PetscPrintf(appctx->comm,"Computed solution vector\n");CHKERRQ(ierr);
ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
ierr = PetscPrintf(appctx->comm,"Exact solution vector\n");CHKERRQ(ierr);
ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
}
/*
Compute the 2-norm and max-norm of the error
*/
ierr = VecAXPY(appctx->solution,-1.0,u);CHKERRQ(ierr);
ierr = VecNorm(appctx->solution,NORM_2,&en2);CHKERRQ(ierr);
en2s = PetscSqrtReal(appctx->h)*en2; /* scale the 2-norm by the grid spacing */
ierr = VecNorm(appctx->solution,NORM_MAX,&enmax);CHKERRQ(ierr);
/*
PetscPrintf() causes only the first processor in this
communicator to print the timestep information.
*/
ierr = PetscPrintf(appctx->comm,"Timestep %D: time = %G,2-norm error = %G, max norm error = %G\n",
step,time,en2s,enmax);CHKERRQ(ierr);
/*
Print debugging information if desired
*/
/* if (appctx->debug) {
ierr = PetscPrintf(appctx->comm,"Error vector\n");CHKERRQ(ierr);
ierr = VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
} */
return 0;
}示例14: test_view
PetscErrorCode test_view( void )
{
Vec X, a,b;
Vec c,d,e,f;
Vec tmp_buf[2];
IS tmp_is[2];
PetscInt index;
PetscReal val;
PetscInt list[]={0,1,2};
PetscScalar vals[]={0.720032,0.061794,0.0100223};
PetscErrorCode ierr;
PetscBool explcit = PETSC_FALSE;
PetscFunctionBegin;
PetscPrintf( PETSC_COMM_WORLD, "\n\n============== %s ==============\n", PETSC_FUNCTION_NAME );
ierr = VecCreate( PETSC_COMM_WORLD, &c );CHKERRQ(ierr);
ierr = VecSetSizes( c, PETSC_DECIDE, 3 );CHKERRQ(ierr);
ierr = VecSetFromOptions( c );CHKERRQ(ierr);
ierr = VecDuplicate( c, &d );CHKERRQ(ierr);
ierr = VecDuplicate( c, &e );CHKERRQ(ierr);
ierr = VecDuplicate( c, &f );CHKERRQ(ierr);
ierr = VecSet( c, 1.0 );CHKERRQ(ierr);
ierr = VecSet( d, 2.0 );CHKERRQ(ierr);
ierr = VecSet( e, 3.0 );CHKERRQ(ierr);
ierr = VecSetValues(f,3,list,vals,INSERT_VALUES);CHKERRQ(ierr);
ierr = VecAssemblyBegin(f);CHKERRQ(ierr);
ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
ierr = VecScale( f, 10.0 );CHKERRQ(ierr);
tmp_buf[0] = e;
tmp_buf[1] = f;
ierr = PetscOptionsGetBool(0,"-explicit_is",&explcit,0);CHKERRQ(ierr);
ierr = GetISs(tmp_buf,tmp_is);CHKERRQ(ierr);
ierr = VecCreateNest(PETSC_COMM_WORLD,2,explcit?tmp_is:PETSC_NULL,tmp_buf,&b);CHKERRQ(ierr);
ierr = VecDestroy(&e);CHKERRQ(ierr);
ierr = VecDestroy(&f);CHKERRQ(ierr);
ierr = ISDestroy(&tmp_is[0]);CHKERRQ(ierr);
ierr = ISDestroy(&tmp_is[1]);CHKERRQ(ierr);
tmp_buf[0] = c;
tmp_buf[1] = d;
ierr = VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&a);CHKERRQ(ierr);
ierr = VecDestroy(&c);CHKERRQ(ierr); ierr = VecDestroy(&d);CHKERRQ(ierr);
tmp_buf[0] = a;
tmp_buf[1] = b;
ierr = VecCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL,tmp_buf,&X);CHKERRQ(ierr);
ierr = VecDestroy(&a);CHKERRQ(ierr);
ierr = VecAssemblyBegin(X);CHKERRQ(ierr);
ierr = VecAssemblyEnd(X);CHKERRQ(ierr);
ierr = VecMax( b, &index, &val );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index );
ierr = VecMin( b, &index, &val );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(min-b) = %f : index = %d \n", val, index );
ierr = VecDestroy(&b);CHKERRQ(ierr);
ierr = VecMax( X, &index, &val );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(max-X) = %f : index = %d \n", val, index );
ierr = VecMin( X, &index, &val );CHKERRQ(ierr);
PetscPrintf( PETSC_COMM_WORLD, "(min-X) = %f : index = %d \n", val, index );
ierr = VecView( X, PETSC_VIEWER_STDOUT_WORLD );CHKERRQ(ierr);
ierr = VecDestroy(&X);CHKERRQ(ierr);
PetscFunctionReturn(0);
}示例15: FormJacobianLocal
/*
FormFunctionLocal - Form the local residual F from the local input X
Input Parameters:
+ dm - The mesh
. X - Local input vector
- user - The user context
Output Parameter:
. F - Local output vector
Note:
We form the residual one batch of elements at a time. This allows us to offload work onto an accelerator,
like a GPU, or vectorize on a multicore machine.
.seealso: FormJacobianLocal()
*/
PetscErrorCode FormFunctionLocal(DM dm, Vec X, Vec F, AppCtx *user)
{
const PetscInt debug = user->debug;
const PetscInt dim = user->dim;
PetscReal *coords, *v0, *J, *invJ, *detJ;
PetscScalar *elemVec, *u;
const PetscInt numCells = cEnd - cStart;
PetscInt cellDof = 0;
PetscInt maxQuad = 0;
PetscInt jacSize = 1;
PetscInt cStart, cEnd, c, field;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = PetscLogEventBegin(user->residualEvent,0,0,0,0);CHKERRQ(ierr);
ierr = VecSet(F, 0.0);CHKERRQ(ierr);
ierr = DMDAGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr);
for(field = 0; field < numFields; ++field) {
PetscInt dof = 1;
for(d = 0; d < dim; ++d) {dof *= user->q[field].numBasisFuncs*user->q[field].numComponents;}
cellDof += dof;
maxQuad = PetscMax(maxQuad, user->q[field].numQuadPoints);
}
for(d = 0; d < dim; ++d) {jacSize *= maxQuad;}
ierr = PetscMalloc3(dim,PetscReal,&coords,dim,PetscReal,&v0,jacSize,PetscReal,&J);CHKERRQ(ierr);
ierr = PetscMalloc4(numCells*cellDof,PetscScalar,&u,numCells*jacSize,PetscReal,&invJ,numCells,PetscReal,&detJ,numCells*cellDof,PetscScalar,&elemVec);CHKERRQ(ierr);
for(c = cStart; c < cEnd; ++c) {
const PetscScalar *x;
PetscInt i;
ierr = DMDAComputeCellGeometry(dm, c, v0, J, &invJ[c*jacSize], &detJ[c]);CHKERRQ(ierr);
if (detJ[c] <= 0.0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Invalid determinant %g for element %d", detJ[c], c);
ierr = DMDAVecGetClosure(dm, PETSC_NULL, X, c, &x);CHKERRQ(ierr);
for(i = 0; i < cellDof; ++i) {
u[c*cellDof+i] = x[i];
}
}
for(field = 0; field < numFields; ++field) {
const PetscInt numQuadPoints = user->q[field].numQuadPoints;
const PetscInt numBasisFuncs = user->q[field].numBasisFuncs;
void (*f0)(PetscScalar u[], const PetscScalar gradU[], PetscScalar f0[]) = user->f0Funcs[field];
void (*f1)(PetscScalar u[], const PetscScalar gradU[], PetscScalar f1[]) = user->f1Funcs[field];
/* Conforming batches */
PetscInt blockSize = numBasisFuncs*numQuadPoints;
PetscInt numBlocks = 1;
PetscInt batchSize = numBlocks * blockSize;
PetscInt numBatches = user->numBatches;
PetscInt numChunks = numCells / (numBatches*batchSize);
ierr = IntegrateResidualBatchCPU(numChunks*numBatches*batchSize, numFields, field, u, invJ, detJ, user->q, f0, f1, elemVec, user);CHKERRQ(ierr);
/* Remainder */
PetscInt numRemainder = numCells % (numBatches * batchSize);
PetscInt offset = numCells - numRemainder;
ierr = IntegrateResidualBatchCPU(numRemainder, numFields, field, &u[offset*cellDof], &invJ[offset*dim*dim], &detJ[offset],
user->q, f0, f1, &elemVec[offset*cellDof], user);CHKERRQ(ierr);
}
for(c = cStart; c < cEnd; ++c) {
if (debug) {ierr = DMPrintCellVector(c, "Residual", cellDof, &elemVec[c*cellDof]);CHKERRQ(ierr);}
ierr = DMComplexVecSetClosure(dm, PETSC_NULL, F, c, &elemVec[c*cellDof], ADD_VALUES);CHKERRQ(ierr);
}
ierr = PetscFree4(u,invJ,detJ,elemVec);CHKERRQ(ierr);
ierr = PetscFree3(coords,v0,J);CHKERRQ(ierr);
if (user->showResidual) {
PetscInt p;
ierr = PetscPrintf(PETSC_COMM_WORLD, "Residual:\n");CHKERRQ(ierr);
for(p = 0; p < user->numProcs; ++p) {
if (p == user->rank) {
Vec f;
ierr = VecDuplicate(F, &f);CHKERRQ(ierr);
ierr = VecCopy(F, f);CHKERRQ(ierr);
ierr = VecChop(f, 1.0e-10);CHKERRQ(ierr);
ierr = VecView(f, PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr);
ierr = VecDestroy(&f);CHKERRQ(ierr);
}
ierr = PetscBarrier((PetscObject) dm);CHKERRQ(ierr);
}
}
ierr = PetscLogEventEnd(user->residualEvent,0,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}