System Preprocessors
u14.c
Go to the documentation of this file.
00001 #include <stdlib.h>
00002 #include "syspro.h"
00003 #include "sysprotransform.h"
00004 #include "sysprolinear.h"
00005 
00006 /* test for the syspro linear package
00007    one preprocessor used: scaling
00008    this uses PETSc's KSP for actually solving a linear system
00009 */
00010 
00011 /* this example features a preprocessor that we can loop over:
00012    preconditioners
00013 */
00014 #undef __FUNCT__
00015 #define __FUNCT__ "create_solver"
00016 /*! Create a solver and 
00017   install a monitor that dynamically increases the maximum
00018   number of iterations.
00019 */
00020 static PetscErrorCode create_solver(NumericalProblem prob,void **ctx)
00021 {
00022   MPI_Comm comm; KSP solver; PetscErrorCode ierr;
00023   PetscFunctionBegin;
00024   ierr = NumericalProblemGetComm(prob,&comm); CHKERRQ(ierr);
00025   ierr = KSPCreate(comm,&solver); CHKERRQ(ierr);
00026   *(KSP*)ctx = solver;
00027   PetscFunctionReturn(0);
00028 }
00029 
00030 #undef __FUNCT__
00031 #define __FUNCT__ "destroy_solver"
00032 static PetscErrorCode destroy_solver(void *ctx)
00033 {
00034   KSP solver; PetscErrorCode ierr;
00035   PetscFunctionBegin;
00036   solver = (KSP)ctx;
00037   ierr = KSPDestroy(&solver); CHKERRQ(ierr);
00038   PetscFunctionReturn(0);
00039 }
00040 
00041 #undef __FUNCT__
00042 #define __FUNCT__ "setup_pc_choices"
00043 static PetscErrorCode setup_pc_choices()
00044 {
00045   PetscErrorCode ierr;
00046   PetscFunctionBegin;
00047   ierr = NewTransformObject("pc",PCNONE,PETSC_NULL); CHKERRQ(ierr);
00048   ierr = NewTransformObject("pc",PCJACOBI,PETSC_NULL); CHKERRQ(ierr);
00049   ierr = NewTransformObject("pc",PCILU,PETSC_NULL); CHKERRQ(ierr);
00050   PetscFunctionReturn(0);
00051 }
00052 
00053 #undef __FUNCT__
00054 #define __FUNCT__ "setup_pc"
00055 static PetscErrorCode setup_pc
00056 (const char *type,int pcv,PetscBool overwrite,
00057  NumericalProblem inproblem,NumericalProblem *outproblem,
00058  void *gctx,void **ctx,PetscBool *success)
00059 {
00060   LinearSystem problem = (LinearSystem)inproblem;
00061   KSP solver; PC pc;  Mat A;
00062   PetscErrorCode ierr;
00063 
00064   PetscFunctionBegin;
00065   *outproblem = inproblem;
00066 
00067   /*
00068    * the fun part is setting up the Petsc PC
00069    */
00070   ierr = LinearSystemGetParts
00071     (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00072   solver = (KSP)gctx;
00073   ierr = KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00074   ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
00075   ierr = PCSetType(pc,type); CHKERRQ(ierr);
00076   *success = PETSC_TRUE;
00077 
00078   PetscFunctionReturn(0);
00079 }
00080 
00081 #undef __FUNCT__
00082 #define __FUNCT__ "unset_pc"
00083 static PetscErrorCode unset_pc
00084 (const char *type,PetscBool overwrite,void *ctx,void *gctx,
00085  NumericalProblem thisproblem,NumericalProblem upproblem,
00086  NumericalSolution old,NumericalSolution nnew)
00087 {
00088 #if 0
00089   LinearSystem problem = (LinearSystem)thisproblem;
00090   Mat A;
00091   PetscErrorCode ierr;
00092 #endif
00093   PetscFunctionBegin;
00094 
00095 #if 0
00096   ierr = LinearSolutionCopy
00097     ((LinearSolution)old,(LinearSolution)nnew); CHKERRQ(ierr);
00098 #endif
00099 
00100   PetscFunctionReturn(0);
00101 }
00102 
00103 
00104 #undef __FUNCT__
00105 #define __FUNCT__ "solvelinear"
00106 static PetscErrorCode solvelinear
00107 (NumericalProblem problem,void *dum,NumericalSolution *rsol)
00108 {
00109   LinearSystem linear = (LinearSystem)problem;
00110   KSP solver; Mat A; Vec rhs,sol; 
00111   KSPConvergedReason reason; PetscInt it; PetscErrorCode ierr;
00112   PetscFunctionBegin;
00113   ierr = LinearSystemGetParts(linear,&A,NULL,&rhs,NULL,NULL); CHKERRQ(ierr);
00114   ierr = VecDuplicate(rhs,&sol); CHKERRQ(ierr);
00115   ierr = KSPCreate(MPI_COMM_SELF,&solver); CHKERRQ(ierr);
00116   ierr = KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
00117   ierr = KSPSolve(solver,rhs,sol); CHKERRQ(ierr);
00118   ierr = KSPGetConvergedReason(solver,&reason); CHKERRQ(ierr);
00119   if (reason<=0) SETERRQ(PETSC_COMM_SELF,1,"This should have converged");
00120   ierr = KSPGetIterationNumber(solver,&it); CHKERRQ(ierr);
00121   printf("Convergence in %d iterations\n",it);
00122   ierr = KSPDestroy(&solver); CHKERRQ(ierr);
00123   ierr = LinearCreateNumericalSolution(NULL,rsol); CHKERRQ(ierr);
00124   ierr = LinearSolutionSetVector((LinearSolution)*rsol,sol); CHKERRQ(ierr);
00125   PetscFunctionReturn(0);
00126 }
00127 
00128 int main(int argc,char **argv) {
00129   NumericalProblem problem; NumericalSolution solution;
00130   Mat A; Vec rhs,sol; PetscErrorCode ierr;
00131   PetscInitialize(&argc,&argv,0,0);
00132   ierr = SysProInitialize(); CHKERRQ(ierr);
00133 
00134   /* set up a simple problem */
00135   {
00136     LinearSystem linearproblem;
00137     int n=100,i; PetscScalar two=2.,mone=-1.,one=1.;
00138 
00139 #include "testmat.c"
00140     /*
00141      * now create the linear system to contain them
00142      */
00143     ierr = CreateLinearSystem(MPI_COMM_SELF,&linearproblem); CHKERRQ(ierr);
00144     ierr = LinearSystemSetParts
00145       (linearproblem,A,NULL,rhs,sol,NULL); CHKERRQ(ierr);
00146     problem = (NumericalProblem)linearproblem;
00147   }
00148 
00149   /* declare the scalings */
00150   ierr = DeclareScalingPreprocessor(); CHKERRQ(ierr);
00151   ierr = PetscOptionsSetValue("-syspro_scaling","exhaustive"); CHKERRQ(ierr);
00152 
00153   /* declare a preconditioner loop */
00154   ierr = DeclarePreprocessor
00155     (/* name */ "pc",
00156      /* one-time creation */ &setup_pc_choices,
00157      /* creation and deletion for specific choices */ NULL,NULL,NULL,
00158      /* create/delete context that will live for the duration of this
00159         preprocessor */
00160      &create_solver,&destroy_solver,
00161      /* set and unset general/specific ? choice */
00162      &setup_pc,NULL /*&unset_pc*/ ); CHKERRQ(ierr);
00163   /* and make sure you loop over it */
00164   ierr = PetscOptionsSetValue("-syspro_pc","none"); CHKERRQ(ierr);
00165 
00166   /* declare functions */
00167   ierr = SysProDeclareFunctions
00168     (NULL,NULL,NULL, /* general setup functions */
00169      solvelinear,LinearDeleteNumericalProblem, /* problem functions */
00170      &LinearCreateNumericalSolution,NULL,LinearDeleteNumericalSolution, /* solution functions */
00171      NULL,NULL,NULL /* context functions */
00172      ); CHKERRQ(ierr);
00173 
00174   ierr = PreprocessorsOptionsHandling(); CHKERRQ(ierr);
00175 
00176   /* give some feedback on what's happening */
00177   ierr = SysProDeclareTraceFunction(SysProDefaultTrace); CHKERRQ(ierr);
00178 
00179   /* solve this problem */
00180   ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00181 
00182   /* check for correctness */
00183   {
00184     Vec sol,asol; int mone=-1; PetscReal nrm;
00185     ierr = LinearSolutionGetVector
00186       ((LinearSolution)solution,&sol); CHKERRQ(ierr);
00187     ierr = VecDuplicate(sol,&asol); CHKERRQ(ierr);
00188     ierr = MatMult(A,sol,asol); CHKERRQ(ierr);
00189     ierr = VecAXPY(asol,mone,rhs); CHKERRQ(ierr);
00190     ierr = VecNorm(asol,NORM_2,&nrm); CHKERRQ(ierr);
00191     if (nrm>1.e-6)
00192       SETERRQ1(PETSC_COMM_SELF,1,"Computed wrong solution: norm residual is %e",nrm);
00193     ierr = VecDestroy(&asol); CHKERRQ(ierr);
00194   }
00195   ierr = LinearDeleteNumericalSolution(solution); CHKERRQ(ierr);
00196   ierr = MatDestroy(&A); CHKERRQ(ierr);
00197   ierr = VecDestroy(&rhs); CHKERRQ(ierr);
00198   ierr = VecDestroy(&sol); CHKERRQ(ierr);
00199   ierr = PetscFree(problem); CHKERRQ(ierr);
00200   ierr = SysProFinalize(); CHKERRQ(ierr);
00201   PetscFinalize();
00202   PetscFunctionReturn(0);
00203 }