00001 #include <stdlib.h>
00002 #include "syspro.h"
00003 #include "anamodsalsamodules.h"
00004 #include "sysprotransform.h"
00005 #include "sysprolinear.h"
00006 #include "anamod.h"
00007 #include "nmd.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017 #undef __FUNCT__
00018 #define __FUNCT__ "create_solver"
00019
00020
00021
00022
00023 static PetscErrorCode create_solver(NumericalProblem prob,void **ctx)
00024 {
00025 MPI_Comm comm; KSP solver; PetscErrorCode ierr;
00026 PetscFunctionBegin;
00027 ierr = NumericalProblemGetComm(prob,&comm); CHKERRQ(ierr);
00028 ierr = KSPCreate(comm,&solver); CHKERRQ(ierr);
00029 *(KSP*)ctx = solver;
00030 PetscFunctionReturn(0);
00031 }
00032
00033 #undef __FUNCT__
00034 #define __FUNCT__ "destroy_solver"
00035 static PetscErrorCode destroy_solver(void *ctx)
00036 {
00037 KSP solver; PetscErrorCode ierr;
00038 PetscFunctionBegin;
00039 solver = (KSP)ctx;
00040 ierr = KSPDestroy(solver); CHKERRQ(ierr);
00041 PetscFunctionReturn(0);
00042 }
00043
00044 #undef __FUNCT__
00045 #define __FUNCT__ "setup_pc_choices"
00046 static PetscErrorCode setup_pc_choices()
00047 {
00048 PetscErrorCode ierr;
00049 PetscFunctionBegin;
00050 ierr = NewTransformObject("pc",PCNONE,PETSC_NULL); CHKERRQ(ierr);
00051 ierr = NewTransformObject("pc",PCJACOBI,PETSC_NULL); CHKERRQ(ierr);
00052 ierr = NewTransformObject("pc",PCILU,PETSC_NULL); CHKERRQ(ierr);
00053 PetscFunctionReturn(0);
00054 }
00055
00056 #undef __FUNCT__
00057 #define __FUNCT__ "setup_pc"
00058 static PetscErrorCode setup_pc
00059 (char *type,int pcv,PetscTruth overwrite,
00060 NumericalProblem inproblem,NumericalProblem *outproblem,
00061 void *gctx,void **ctx,PetscTruth *success)
00062 {
00063 LinearSystem problem = (LinearSystem)inproblem;
00064 KSP solver; PC pc; Mat A;
00065 PetscErrorCode ierr;
00066
00067 PetscFunctionBegin;
00068 *outproblem = inproblem;
00069
00070
00071
00072
00073 ierr = LinearSystemGetParts
00074 (problem,&A,PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL); CHKERRQ(ierr);
00075 solver = (KSP)gctx;
00076 ierr = KSPSetOperators(solver,A,A,SAME_NONZERO_PATTERN); CHKERRQ(ierr);
00077 ierr = KSPGetPC(solver,&pc); CHKERRQ(ierr);
00078 ierr = PCSetType(pc,type); CHKERRQ(ierr);
00079 *success = PETSC_TRUE;
00080
00081 PetscFunctionReturn(0);
00082 }
00083
00084 #undef __FUNCT__
00085 #define __FUNCT__ "unset_pc"
00086 static PetscErrorCode unset_pc
00087 (char *type,PetscTruth overwrite,void *ctx,void *gctx,
00088 NumericalProblem thisproblem,NumericalProblem upproblem,
00089 NumericalSolution old,NumericalSolution new)
00090 {
00091 #if 0
00092 LinearSystem problem = (LinearSystem)thisproblem;
00093 Mat A;
00094 PetscErrorCode ierr;
00095 #endif
00096 PetscFunctionBegin;
00097
00098 #if 0
00099 ierr = LinearSolutionCopy
00100 ((LinearSolution)old,(LinearSolution)new); CHKERRQ(ierr);
00101 #endif
00102
00103 PetscFunctionReturn(0);
00104 }
00105
00106
00107 #undef __FUNCT__
00108 #define __FUNCT__ "solvelinear"
00109 static PetscErrorCode solvelinear
00110 (NumericalProblem problem,void *dum,NumericalSolution *rsol)
00111 {
00112 LinearSystem linear = (LinearSystem)problem;
00113 KSP solver; Mat A; Vec rhs,sol;
00114 KSPConvergedReason reason; PetscInt it; PetscErrorCode ierr;
00115 PetscFunctionBegin;
00116 ierr = LinearSystemGetParts(linear,&A,NULL,&rhs,NULL,NULL); CHKERRQ(ierr);
00117 ierr = VecDuplicate(rhs,&sol); CHKERRQ(ierr);
00118 ierr = KSPCreate(MPI_COMM_SELF,&solver); CHKERRQ(ierr);
00119 ierr = KSPSetOperators(solver,A,A,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);
00120 ierr = KSPSolve(solver,rhs,sol); CHKERRQ(ierr);
00121 ierr = KSPGetConvergedReason(solver,&reason); CHKERRQ(ierr);
00122 if (reason<=0) SETERRQ(1,"This should have converged");
00123 ierr = KSPGetIterationNumber(solver,&it); CHKERRQ(ierr);
00124 printf("Convergence in %d iterations\n",it);
00125 ierr = KSPDestroy(solver); CHKERRQ(ierr);
00126 ierr = LinearCreateNumericalSolution(NULL,rsol); CHKERRQ(ierr);
00127 ierr = LinearSolutionSetVector((LinearSolution)*rsol,sol); CHKERRQ(ierr);
00128 PetscFunctionReturn(0);
00129 }
00130
00131 int main(int argc,char **argv) {
00132 NumericalProblem problem; NumericalSolution solution;
00133 Mat A; Vec rhs,sol; PetscErrorCode ierr;
00134 PetscInitialize(&argc,&argv,0,0);
00135 ierr = SysProInitialize(); CHKERRQ(ierr);
00136
00137
00138 {
00139 LinearSystem linearproblem; NMD_metadata nmd;
00140 int n=100,i; PetscScalar two=2.,mone=-1.,one=1.;
00141
00142 #include "testmat16.c"
00143
00144
00145
00146 ierr = CreateLinearSystem(MPI_COMM_SELF,&linearproblem); CHKERRQ(ierr);
00147 ierr = NMDCreateObject(&nmd); CHKERRQ(ierr);
00148 ierr = RegisterStructureModules(); CHKERRQ(ierr);
00149 ierr = LinearSystemSetMetadata(linearproblem,nmd); CHKERRQ(ierr);
00150 ierr = LinearSystemSetParts
00151 (linearproblem,A,NULL,rhs,sol,NULL); CHKERRQ(ierr);
00152 problem = (NumericalProblem)linearproblem;
00153 }
00154
00155
00156 ierr = DeclareSingletonPreprocessor(); CHKERRQ(ierr);
00157
00158
00159 ierr = DeclarePreprocessor
00160 ( "pc",
00161 &setup_pc_choices,
00162 NULL,NULL,NULL,
00163
00164
00165 &create_solver,&destroy_solver,
00166
00167 &setup_pc,NULL ); CHKERRQ(ierr);
00168
00169 ierr = PetscOptionsSetValue("-syspro_pc","jacobi"); CHKERRQ(ierr);
00170
00171
00172 ierr = SysProDeclareFunctions
00173 (NULL,NULL,NULL,
00174 solvelinear,NULL,
00175 &LinearCreateNumericalSolution,NULL,LinearDeleteNumericalSolution,
00176 NULL,NULL,NULL
00177 ); CHKERRQ(ierr);
00178
00179 ierr = PreprocessorsOptionsHandling(); CHKERRQ(ierr);
00180
00181
00182 ierr = SysProDeclareTraceFunction(SysProDefaultTrace); CHKERRQ(ierr);
00183
00184
00185 ierr = PreprocessedProblemSolving(problem,&solution); CHKERRQ(ierr);
00186
00187
00188 {
00189 Vec sol,asol; int mone=-1; PetscReal nrm;
00190 ierr = LinearSolutionGetVector
00191 ((LinearSolution)solution,&sol); CHKERRQ(ierr);
00192 ierr = VecDuplicate(sol,&asol); CHKERRQ(ierr);
00193 ierr = MatMult(A,sol,asol); CHKERRQ(ierr);
00194 ierr = VecAXPY(asol,mone,rhs); CHKERRQ(ierr);
00195 ierr = VecNorm(asol,NORM_2,&nrm); CHKERRQ(ierr);
00196 if (nrm>1.e-6)
00197 SETERRQ1(1,"Computed wrong solution: norm residual is %e",nrm);
00198 ierr = VecDestroy(asol); CHKERRQ(ierr);
00199 }
00200
00201 {
00202 NMD_metadata nmd;
00203 ierr = LinearSystemGetMetadata((LinearSystem)problem,&nmd); CHKERRQ(ierr);
00204 ierr = NMDDestroyObject(nmd); CHKERRQ(ierr);
00205 ierr = DeregisterModules(); CHKERRQ(ierr);
00206 }
00207
00208 ierr = LinearDeleteNumericalSolution(solution); CHKERRQ(ierr);
00209 ierr = MatDestroy(A); CHKERRQ(ierr);
00210 ierr = VecDestroy(rhs); CHKERRQ(ierr);
00211 ierr = VecDestroy(sol); CHKERRQ(ierr);
00212 ierr = PetscFree(problem); CHKERRQ(ierr);
00213 ierr = SysProFinalize(); CHKERRQ(ierr);
00214 PetscFinalize();
00215 PetscFunctionReturn(0);
00216 }