Attached is the source code and the makefile used to compile/run the code. The source code is basically a dumbed down version of SNES ex12 plus necessary TAO routines.
Thanks On Tue, Oct 14, 2014 at 8:39 PM, Matthew Knepley <knep...@gmail.com> wrote: > On Tue, Oct 14, 2014 at 8:04 PM, Justin Chang <jychan...@gmail.com> wrote: > >> Hi all, >> >> So I am writing a non-negative diffusion equation using DMPlex's FEM and >> Tao's SetVariableBounds functions. My code works really perfectly when I >> run it with one processor. However, once I use 2 or more processors, I get >> this error: >> > > It looks like the problem is in the TAO definition, but you can check by > just solving your problem with, for instance, BJacobi-LU in parallel. > > >> [0]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> [0]PETSC ERROR: Nonconforming object sizes >> [0]PETSC ERROR: Vector wrong size 89 for scatter 88 (scatter reverse and >> vector to != ctx from size) >> [1]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> [1]PETSC ERROR: Nonconforming object sizes >> [1]PETSC ERROR: Vector wrong size 87 for scatter 88 (scatter reverse and >> vector to != ctx from size) >> [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >> for trouble shooting. >> [1]PETSC ERROR: Petsc Development GIT revision: v3.5.2-526-gfaecc80 GIT >> Date: 2014-10-04 20:10:35 -0500 >> [1]PETSC ERROR: [0]PETSC ERROR: See >> http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. >> [0]PETSC ERROR: Petsc Development GIT revision: v3.5.2-526-gfaecc80 GIT >> Date: 2014-10-04 20:10:35 -0500 >> [0]PETSC ERROR: ./bin/diff2D on a arch-linux2-c-debug named pacotaco by >> justin Tue Oct 14 19:48:50 2014 >> [0]PETSC ERROR: ./bin/diff2D on a arch-linux2-c-debug named pacotaco by >> justin Tue Oct 14 19:48:50 2014 >> [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >> --with-fc=gfortran --download-fblaslapack --download-mpich >> --with-debugging=1 --download-metis --download-parmetis --download-triangle >> --with-cmake=cmake --download-ctetgen --download-superlu >> --download-scalapack --download-mumps --download-hdf5 --with-valgrind=1 >> -with-cmake=cmake >> [1]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >> --with-fc=gfortran --download-fblaslapack --download-mpich >> --with-debugging=1 --download-metis --download-parmetis --download-triangle >> --with-cmake=cmake --download-ctetgen --download-superlu >> --download-scalapack --download-mumps --download-hdf5 --with-valgrind=1 >> -with-cmake=cmake >> [0]PETSC ERROR: #1 VecScatterBegin() line 1713 in >> /home/justin/petsc-master/src/vec/vec/utils/vscat.c >> #1 VecScatterBegin() line 1713 in >> /home/justin/petsc-master/src/vec/vec/utils/vscat.c >> [1]PETSC ERROR: [0]PETSC ERROR: #2 MatMultTranspose_MPIAIJ() line 1010 in >> /home/justin/petsc-master/src/mat/impls/aij/mpi/mpiaij.c >> [0]PETSC ERROR: #2 MatMultTranspose_MPIAIJ() line 1010 in >> /home/justin/petsc-master/src/mat/impls/aij/mpi/mpiaij.c >> [1]PETSC ERROR: #3 MatMultTranspose() line 2242 in >> /home/justin/petsc-master/src/mat/interface/matrix.c >> #3 MatMultTranspose() line 2242 in >> /home/justin/petsc-master/src/mat/interface/matrix.c >> [0]PETSC ERROR: #4 IPMComputeKKT() line 616 in >> /home/justin/petsc-master/src/tao/constrained/impls/ipm/ipm.c >> [1]PETSC ERROR: #4 IPMComputeKKT() line 616 in >> /home/justin/petsc-master/src/tao/constrained/impls/ipm/ipm.c >> [1]PETSC ERROR: [0]PETSC ERROR: #5 TaoSolve_IPM() line 50 in >> /home/justin/petsc-master/src/tao/constrained/impls/ipm/ipm.c >> [0]PETSC ERROR: #5 TaoSolve_IPM() line 50 in >> /home/justin/petsc-master/src/tao/constrained/impls/ipm/ipm.c >> [1]PETSC ERROR: #6 TaoSolve() line 190 in >> /home/justin/petsc-master/src/tao/interface/taosolver.c >> #6 TaoSolve() line 190 in >> /home/justin/petsc-master/src/tao/interface/taosolver.c >> [0]PETSC ERROR: #7 main() line 341 in >> /home/justin/Dropbox/Research_Topics/Petsc_Nonneg_diffusion/src/diff2D.c >> [1]PETSC ERROR: #7 main() line 341 in >> /home/justin/Dropbox/Research_Topics/Petsc_Nonneg_diffusion/src/diff2D.c >> [1]PETSC ERROR: [0]PETSC ERROR: ----------------End of Error Message >> -------send entire error message to petsc-ma...@mcs.anl.gov---------- >> ----------------End of Error Message -------send entire error message to >> petsc-ma...@mcs.anl.gov---------- >> application called MPI_Abort(MPI_COMM_WORLD, 60) - process 0 >> application called MPI_Abort(MPI_COMM_WORLD, 60) - process 1 >> [cli_1]: aborting job: >> application called MPI_Abort(MPI_COMM_WORLD, 60) - process 1 >> [cli_0]: aborting job: >> application called MPI_Abort(MPI_COMM_WORLD, 60) - process 0 >> >> >> I have no idea how or why I am getting this error. What does this mean? >> > > It looks like one dof is given to proc 0 which should live on proc 1. We > have to look at the divisions > made in the KKT solver. Can you send a small example? > > Matt > > >> My code is essentially built off of SNES ex12.c. The Jacobian matrix, >> residual vector, and solution vector were created using DMPlex and the >> built-in FEM functions within. The Hessian matrix and gradient vector were >> created by simple matmult() functions of the jacobian and residual. The >> lower bounds vector was created by duplicating the solution vector (initial >> guess set to zero). My FormFunctionGradient() is basically the same thing >> as in the maros.c example. can give more information if needed. >> >> Thanks, >> Justin >> > > > > -- > What most experimenters take for granted before they begin their > experiments is infinitely more interesting than any results to which their > experiments lead. > -- Norbert Wiener >
static char help[] = "Nonnegative diffusion 2D using P1 triangular elements";
#include <petscdmplex.h>
#include <petscsnes.h>
#include <petscds.h>
//#include <petscviewerhdf5.h>
#include <petscmath.h>
#include <petsctao.h>
PetscInt spatialDim = 2;
typedef struct {
PetscLogEvent createMeshEvent;
char meshfile[2048]; /* The optional mesh file */
PetscReal refinementLimit; /* The largest allowable cell volume */
PetscBool flg_exo,flg_gmsh; /* Flags for exodus or GMSH mesh */
PetscBool nonneg; /* Flag for solving nonnegative equation */
PetscViewer viewer; /* Multi-purpose viewer */
Mat A,J; /* Jacobian/Hessian matrices */
Vec u,r,b; /* Solution and residual vectors */
Vec f,lb; /* Gradient and lower bounds vectors */
void (**bcFuncs)(const PetscReal x[], PetscScalar *u, void *ctx);
} AppCtx;
/* Initial guess */
void zero(const PetscReal coords[], PetscScalar *u, void *ctx)
{
*u = 0.0;
}
/* Boundary conditions */
void bcs(const PetscReal x[], PetscScalar *u, void *ctx)
{
*u = 0.0;
}
/* Residual -<w,f(x)> */
void f0_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[])
{
if (x[0] >= 3.0/8.0 && x[0] <= 5.0/8.0 && x[1] >= 3.0/8.0 && x[1] <= 5.0/8.0)
f0[0] = -1;
else
f0[0] = 0;
}
/* Residual <grad[w],D(x)grad[u]> */
void f1_u(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[])
{
PetscInt d,d2;
for (d = 0; d < spatialDim; d++) {
f1[d] = 0.0;
for (d2 = 0; d2 < spatialDim; d2++)
f1[d] += a[d*spatialDim+d2]*u_x[d2];
}
}
/* Jacobian <grad[w],D(x)grad[u]> */
void g3_uu(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[])
{
PetscInt d,d2;
for (d = 0; d < spatialDim; d++) {
for (d2 = 0; d2 < spatialDim; d2++)
g3[d*spatialDim+d2] = a[d*spatialDim+d2];
}
}
/* Diffusivity tensor */
void aux(const PetscReal x[], PetscScalar *u, void *ctx)
{
PetscScalar eps;
eps = 0.05;
u[0] = x[1]*x[1]+eps*x[0]*x[0];
u[1] = -(1-eps)*x[0]*x[1];
u[2] = -(1-eps)*x[0]*x[1];
u[3] = x[0]*x[0]+eps*x[1]*x[1];
}
#undef __FUNCT__
#define __FUNCT__ "ProcessOptions"
PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
{
PetscErrorCode ierr;
PetscFunctionBeginUser;
options->meshfile[0] = '\0';
options->refinementLimit = 0.0125;
options->flg_exo = PETSC_FALSE;
options->flg_gmsh = PETSC_FALSE;
options->nonneg = PETSC_FALSE;
ierr = PetscOptionsBegin(comm, "", "Diffusion Problem Options", "DMPLEX");CHKERRQ(ierr);
/* Nonnegative optimization */
ierr = PetscOptionsBool("-nonnegative", "Solve with the non-negative constraint", "ex12.c", options->nonneg, &options->nonneg, NULL);CHKERRQ(ierr);
/* Refinement for box mesh only */
ierr = PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "diff2D.c", options->refinementLimit, &options->refinementLimit, NULL);CHKERRQ(ierr);
/* Optional: read exodusii mesh file */
ierr = PetscOptionsString("-exodus", "Exodus.II filename to read", "diff2D.c", options->meshfile, options->meshfile, sizeof(options->meshfile), &options->flg_exo);CHKERRQ(ierr);
#if !defined(PETSC_HAVE_EXODUSII)
if (options->flg_exo) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "This option requires ExodusII support. Reconfigure using --download-exodusii");
#endif
/* Optional: read exodusii mesh file */
ierr = PetscOptionsString("-gmsh", "GMSH filename to read", "diff2D.c", options->meshfile, options->meshfile, sizeof(options->meshfile), &options->flg_gmsh);CHKERRQ(ierr);
if (options->flg_exo && options->flg_gmsh) SETERRQ(comm,PETSC_ERR_ARG_IDN, "Cannot read both Exodusii and GMSH file");
ierr = PetscOptionsEnd();
ierr = PetscLogEventRegister("CreateMesh", DM_CLASSID, &options->createMeshEvent);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "CreateMesh"
PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
{
PetscReal refinementLimit = user->refinementLimit;
const char *meshfile = user->meshfile;
size_t len;
PetscErrorCode ierr;
PetscFunctionBeginUser;
ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
ierr = PetscStrlen(meshfile, &len);CHKERRQ(ierr);
/* Box mesh */
if (!len) {
ierr = DMPlexCreateBoxMesh(comm, spatialDim, PETSC_TRUE, dm);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) *dm, "Mesh");CHKERRQ(ierr);
/* Exodusii file */
} else if (user->flg_exo) {
ierr = DMPlexCreateExodusFromFile(comm, meshfile, PETSC_FALSE, dm);CHKERRQ(ierr);
ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr);
/* GMSH file */
} else if (user->flg_gmsh) {
DMLabel label;
ierr = PetscViewerCreate(comm,&user->viewer);CHKERRQ(ierr);
ierr = PetscViewerSetType(user->viewer,PETSCVIEWERASCII);CHKERRQ(ierr);
ierr = PetscViewerFileSetMode(user->viewer,FILE_MODE_READ);CHKERRQ(ierr);
ierr = PetscViewerFileSetName(user->viewer,user->meshfile);CHKERRQ(ierr);
ierr = DMPlexCreateGmsh(comm,user->viewer,PETSC_TRUE,dm);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&user->viewer);CHKERRQ(ierr);
ierr = DMPlexGetLabel(*dm,"Face Sets",&label);
if (label) {ierr = DMPlexLabelComplete(*dm,label);CHKERRQ(ierr);}
}
{
DM refinedMesh = NULL;
DM distributedMesh = NULL;
/* Refine mesh using a volume constraint */
ierr = DMPlexSetRefinementLimit(*dm, refinementLimit);CHKERRQ(ierr);
ierr = DMRefine(*dm, comm, &refinedMesh);CHKERRQ(ierr);
if (refinedMesh) {
const char *name;
ierr = PetscObjectGetName((PetscObject) *dm, &name);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) refinedMesh, name);CHKERRQ(ierr);
ierr = DMDestroy(dm);CHKERRQ(ierr);
*dm = refinedMesh;
}
/* Distribute mesh over processes */
ierr = DMPlexDistribute(*dm, "metis", 0, NULL, &distributedMesh);CHKERRQ(ierr);
if (distributedMesh) {
ierr = DMDestroy(dm);CHKERRQ(ierr);
*dm = distributedMesh;
}
}
ierr = DMSetFromOptions(*dm);CHKERRQ(ierr);
ierr = DMViewFromOptions(*dm, NULL, "-dm_view");CHKERRQ(ierr);
ierr = PetscLogEventEnd(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "SetupProblem"
PetscErrorCode SetupProblem(DM dm, AppCtx *user)
{
DM cdm = dm;
const PetscInt id = user->flg_gmsh ? 0 : 1;
PetscFE feAux;
PetscFE fe;
PetscDS prob;
PetscErrorCode ierr;
PetscFunctionBeginUser;
/* Create finite element for diffusion */
ierr = PetscFECreateDefault(dm, spatialDim, 1, PETSC_TRUE, NULL, -1, &fe);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) fe, "diffusion");CHKERRQ(ierr);
PetscQuadrature q;
/* Create finite element for diffusivity tensor */
ierr = PetscFECreateDefault(dm, spatialDim, spatialDim*spatialDim, PETSC_TRUE, "mat_", -1, &feAux);CHKERRQ(ierr);
ierr = PetscFEGetQuadrature(fe, &q);CHKERRQ(ierr);
ierr = PetscFESetQuadrature(feAux, q);CHKERRQ(ierr);
/* Set discretization and boundary conditions for each mesh */
while (cdm) {
ierr = DMGetDS(cdm, &prob);CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe);CHKERRQ(ierr);
/* Set discretization for diffusivity */
DM dmAux;
PetscDS probAux;
Vec nu;
void (*matFuncs[1])(const PetscReal x[], PetscScalar *u, void *ctx) = {aux};
ierr = DMClone(cdm, &dmAux);CHKERRQ(ierr);
ierr = DMPlexCopyCoordinates(cdm, dmAux);CHKERRQ(ierr);
ierr = DMGetDS(dmAux, &probAux);CHKERRQ(ierr);
ierr = PetscDSSetDiscretization(probAux, 0, (PetscObject) feAux);CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject) dm, "dmAux", (PetscObject) dmAux);CHKERRQ(ierr);
ierr = DMCreateLocalVector(dmAux, &nu);CHKERRQ(ierr);
ierr = DMPlexProjectFunctionLocal(dmAux, matFuncs, NULL, INSERT_ALL_VALUES, nu);CHKERRQ(ierr);
ierr = PetscObjectCompose((PetscObject) cdm, "A", (PetscObject) nu);CHKERRQ(ierr);
ierr = VecDestroy(&nu);CHKERRQ(ierr);
ierr = DMDestroy(&dmAux);CHKERRQ(ierr);
/* Set discretization for diffusion */
ierr = PetscDSSetResidual(prob, 0, f0_u, f1_u);CHKERRQ(ierr);
ierr = PetscDSSetJacobian(prob, 0, 0, NULL, NULL, NULL, g3_uu);CHKERRQ(ierr);
user->bcFuncs[0] = bcs;
/* Set boundary conditions */
ierr = DMPlexAddBoundary(cdm, PETSC_TRUE, "wall",user->flg_gmsh ? "Face Sets" : "marker", 0, user->bcFuncs[0], 1, &id, user);CHKERRQ(ierr);
ierr = DMPlexGetCoarseDM(cdm, &cdm);CHKERRQ(ierr);
}
ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
ierr = PetscFEDestroy(&feAux);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "FormFunctionGradient"
PetscErrorCode FormFunctionGradient(Tao tao, Vec x, PetscReal *f, Vec g, void *ctx)
{
AppCtx *user = (AppCtx*)ctx;
PetscScalar xtHx;
PetscErrorCode ierr;
PetscFunctionBegin;
ierr = MatMult(user->A,x,g);CHKERRQ(ierr);
ierr = VecDot(x,g,&xtHx);CHKERRQ(ierr);
ierr = VecDot(x,user->f,f);CHKERRQ(ierr);
*f += 0.5*xtHx;
ierr = VecAXPY(g,1.0,user->f);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "FormHessian"
PetscErrorCode FormHessian(Tao tao, Vec x, Mat H, Mat Hpre, void *ctx)
{
PetscFunctionBegin;
PetscFunctionReturn(0);
}
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc, char **argv)
{
DM dm; /* Problem specification */
SNES snes; /* nonlinear solver */
Tao tao; /* Optimization solver */
KSP ksp; /* linear solver */
TaoConvergedReason reason; /* Convergence reason */
AppCtx user; /* user-defined work context */
PetscInt its; /* iterations for convergence */
PetscErrorCode ierr;
ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
ierr = SNESCreate(PETSC_COMM_WORLD, &snes);CHKERRQ(ierr);
ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
ierr = TaoSetType(tao,TAOIPM);CHKERRQ(ierr);
ierr = CreateMesh(PETSC_COMM_WORLD, &user, &dm);CHKERRQ(ierr);
ierr = SNESSetDM(snes, dm);CHKERRQ(ierr);
ierr = DMSetApplicationContext(dm, &user);CHKERRQ(ierr);
ierr = PetscMalloc(1 * sizeof(void (*)(const PetscReal[], PetscScalar *, void *)), &user.bcFuncs);CHKERRQ(ierr);
ierr = SetupProblem(dm, &user);CHKERRQ(ierr);
/* Create global matrices and vectors */
ierr = DMCreateGlobalVector(dm, &user.u);CHKERRQ(ierr);
ierr = PetscObjectSetName((PetscObject) user.u, "fields");CHKERRQ(ierr);
ierr = VecDuplicate(user.u, &user.r);CHKERRQ(ierr);
ierr = VecDuplicate(user.u, &user.b);CHKERRQ(ierr);
ierr = VecDuplicate(user.u, &user.lb);CHKERRQ(ierr);
ierr = VecDuplicate(user.u, &user.f);CHKERRQ(ierr);
ierr = DMSetMatType(dm,MATAIJ);CHKERRQ(ierr);
ierr = DMCreateMatrix(dm, &user.J);CHKERRQ(ierr);
user.A = user.J;
/* Set residual and jacobian evaluation functions */
ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*)) DMPlexSNESComputeResidualFEM, &user);CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM,Vec,Mat,Mat,void*)) DMPlexSNESComputeJacobianFEM, &user);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes, user.A, user.J, NULL, NULL);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
/* Initial guess and lower bound */
void (*initialGuess[1])(const PetscReal x[], PetscScalar *, void *ctx) = {zero};
ierr = DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_VALUES, user.u);CHKERRQ(ierr);
ierr = DMPlexProjectFunction(dm, initialGuess, NULL, INSERT_ALL_VALUES, user.lb);CHKERRQ(ierr);
/* Solve - original Galerkin */
if (!user.nonneg) {
ierr = PetscPrintf(PETSC_COMM_WORLD, "\nClassical Galerkin formalism:\n");CHKERRQ(ierr);
ierr = SNESSolve(snes, NULL, user.u);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes, &its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD, "Number of SNES iterations = %D\n", its);CHKERRQ(ierr);
/* Solve - nonnegative optimization */
} else {
ierr = PetscPrintf(PETSC_COMM_WORLD, "\nNon-negative optimization:\n");CHKERRQ(ierr);
/* Compute residual and jacobian */
ierr = SNESComputeFunction(snes,user.u,user.r);CHKERRQ(ierr);
ierr = SNESComputeJacobian(snes,user.u,user.A,user.A);CHKERRQ(ierr);
/* Form vector necessary for objective and gradient routine */
ierr = MatMult(user.A,user.u,user.b);CHKERRQ(ierr);
ierr = VecChop(user.r,1.0e-10);CHKERRQ(ierr);
ierr = VecChop(user.b,1.0e-10);CHKERRQ(ierr);
ierr = VecWAXPY(user.f,-1.0,user.b,user.r);CHKERRQ(ierr);
/* Set up Tao solver */
ierr = TaoSetInitialVector(tao,user.u);CHKERRQ(ierr);
ierr = TaoSetVariableBounds(tao,user.lb,NULL);CHKERRQ(ierr);
ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user);CHKERRQ(ierr);
ierr = TaoSetHessianRoutine(tao,user.A,user.A,FormHessian,(void*)&user);CHKERRQ(ierr);
ierr = TaoSetTolerances(tao,1e-9,0,0,0,0);CHKERRQ(ierr);
ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* Solve */
ierr = TaoSolve(tao);CHKERRQ(ierr);
ierr = TaoGetConvergedReason(tao,&reason);CHKERRQ(ierr);
if (reason < 0) {
ierr = PetscPrintf(MPI_COMM_WORLD, "TAO failed to converge.\n");CHKERRQ(ierr);
} else {
ierr = PetscPrintf(MPI_COMM_WORLD, "Optimization terminated with status %D.\n", reason);CHKERRQ(ierr);
}
}
/* Output to VTK *
ierr = PetscViewerCreate(PetscObjectComm((PetscObject)dm),&user.viewer);CHKERRQ(ierr);
ierr = PetscViewerSetFormat(user.viewer,PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr);
ierr = PetscViewerSetType(user.viewer, PETSCVIEWERVTK);CHKERRQ(ierr);
if (user.nonneg) {
ierr = PetscViewerFileSetName(user.viewer, "nonneg.vtk");CHKERRQ(ierr);
} else {
ierr = PetscViewerFileSetName(user.viewer, "orig.vtk");CHKERRQ(ierr);
}
ierr = VecView(user.u, user.viewer);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&user.viewer);CHKERRQ(ierr);
*/
/* Clean up */
if (user.A != user.J) {ierr = MatDestroy(&user.A);CHKERRQ(ierr);}
ierr = TaoDestroy(&tao);CHKERRQ(ierr);
ierr = MatDestroy(&user.J);CHKERRQ(ierr);
ierr = VecDestroy(&user.u);CHKERRQ(ierr);
ierr = VecDestroy(&user.r);CHKERRQ(ierr);
ierr = VecDestroy(&user.b);CHKERRQ(ierr);
ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = DMDestroy(&dm);CHKERRQ(ierr);
ierr = PetscFree(user.bcFuncs);CHKERRQ(ierr);
ierr = PetscFinalize();
return 0;
}
makefile
Description: Binary data
