System Preprocessors
|
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 }