1 #include <petsc/private/pcmgimpl.h> /*I "petscksp.h" I*/ 2 #include <petscdm.h> 3 4 static PetscErrorCode PCMGCreateCoarseSpaceDefault_Private(PC pc, PetscInt level, PCMGCoarseSpaceType cstype, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 5 { 6 PetscBool poly = cstype == PCMG_POLYNOMIAL ? PETSC_TRUE : PETSC_FALSE; 7 PetscErrorCode (**funcs)(PetscInt,PetscReal,const PetscReal[],PetscInt,PetscScalar*,void*); 8 void **ctxs; 9 PetscInt dim, d, Nf, f, k; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 ierr = DMGetCoordinateDim(dm, &dim);CHKERRQ(ierr); 14 ierr = DMGetNumFields(dm, &Nf);CHKERRQ(ierr); 15 if (Nc % dim) SETERRQ2(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "The number of coarse vectors %D must be divisible by the dimension %D", Nc, dim); 16 ierr = PetscMalloc2(Nf, &funcs, Nf, &ctxs);CHKERRQ(ierr); 17 if (!*coarseSpace) {ierr = PetscCalloc1(Nc, coarseSpace);CHKERRQ(ierr);} 18 for (k = 0; k < Nc/dim; ++k) { 19 for (f = 0; f < Nf; ++f) {ctxs[f] = &k;} 20 for (d = 0; d < dim; ++d) { 21 if (!(*coarseSpace)[k*dim+d]) {ierr = DMCreateGlobalVector(dm, &(*coarseSpace)[k*dim+d]);CHKERRQ(ierr);} 22 ierr = DMSetBasisFunction_Internal(Nf, poly, d, funcs);CHKERRQ(ierr); 23 ierr = DMProjectFunction(dm, 0.0, funcs, ctxs, INSERT_ALL_VALUES, (*coarseSpace)[k*dim+d]);CHKERRQ(ierr); 24 } 25 } 26 ierr = PetscFree2(funcs, ctxs);CHKERRQ(ierr); 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode PCMGCreateCoarseSpace_Polynomial(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 31 { 32 PetscErrorCode ierr; 33 34 PetscFunctionBegin; 35 ierr = PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_POLYNOMIAL, dm, ksp, Nc, initialGuess, coarseSpace);CHKERRQ(ierr); 36 PetscFunctionReturn(0); 37 } 38 39 PetscErrorCode PCMGCreateCoarseSpace_Harmonic(PC pc, PetscInt level, DM dm, KSP ksp, PetscInt Nc, const Vec initialGuess[], Vec **coarseSpace) 40 { 41 PetscErrorCode ierr; 42 43 PetscFunctionBegin; 44 ierr = PCMGCreateCoarseSpaceDefault_Private(pc, level, PCMG_HARMONIC, dm, ksp, Nc, initialGuess, coarseSpace);CHKERRQ(ierr); 45 PetscFunctionReturn(0); 46 } 47 48 /* 49 PCMGComputeCoarseSpace_Internal - Compute vectors on level l that must be accurately interpolated. 50 51 Input Parameters: 52 + pc - The PCMG 53 . l - The level 54 . Nc - The size of the space (number of vectors) 55 - cspace - The space from level l-1, or NULL 56 57 Output Parameter: 58 . space - The space which must be accurately interpolated. 59 60 Level: developer 61 62 Note: This space is normally used to adapt the interpolator. 63 64 .seealso: PCMGAdaptInterpolator_Private() 65 */ 66 PetscErrorCode PCMGComputeCoarseSpace_Internal(PC pc, PetscInt l, PCMGCoarseSpaceType cstype, PetscInt Nc, const Vec cspace[], Vec *space[]) 67 { 68 PetscErrorCode (*coarseConstructor)(PC, PetscInt, DM, KSP, PetscInt, const Vec[], Vec*[]); 69 DM dm; 70 KSP smooth; 71 PetscErrorCode ierr; 72 73 PetscFunctionBegin; 74 switch (cstype) { 75 case PCMG_POLYNOMIAL: coarseConstructor = &PCMGCreateCoarseSpace_Polynomial;break; 76 case PCMG_HARMONIC: coarseConstructor = &PCMGCreateCoarseSpace_Harmonic;break; 77 case PCMG_EIGENVECTOR: 78 if (l > 0) {ierr = PCMGGetCoarseSpaceConstructor("BAMG_MEV", &coarseConstructor);CHKERRQ(ierr);} 79 else {ierr = PCMGGetCoarseSpaceConstructor("BAMG_EV", &coarseConstructor);CHKERRQ(ierr);} 80 break; 81 case PCMG_GENERALIZED_EIGENVECTOR: 82 if (l > 0) {ierr = PCMGGetCoarseSpaceConstructor("BAMG_MGEV", &coarseConstructor);CHKERRQ(ierr);} 83 else {ierr = PCMGGetCoarseSpaceConstructor("BAMG_GEV", &coarseConstructor);CHKERRQ(ierr);} 84 break; 85 default: SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle coarse space type %D", cstype); 86 } 87 ierr = PCMGGetSmoother(pc, l, &smooth);CHKERRQ(ierr); 88 ierr = KSPGetDM(smooth, &dm);CHKERRQ(ierr); 89 ierr = (*coarseConstructor)(pc, l, dm, smooth, Nc, cspace, space);CHKERRQ(ierr); 90 PetscFunctionReturn(0); 91 } 92 93 /* 94 PCMGAdaptInterpolator_Internal - Adapt interpolator from level l-1 to level 1 95 96 Input Parameters: 97 + pc - The PCMG 98 . l - The level l 99 . csmooth - The (coarse) smoother for level l-1 100 . fsmooth - The (fine) smoother for level l 101 . Nc - The size of the subspace used for adaptation 102 . cspace - The (coarse) vectors in the subspace for level l-1 103 - fspace - The (fine) vectors in the subspace for level l 104 105 Level: developer 106 107 Note: This routine resets the interpolation and restriction for level l. 108 109 .seealso: PCMGComputeCoarseSpace_Private() 110 */ 111 PetscErrorCode PCMGAdaptInterpolator_Internal(PC pc, PetscInt l, KSP csmooth, KSP fsmooth, PetscInt Nc, Vec cspace[], Vec fspace[]) 112 { 113 PC_MG *mg = (PC_MG *) pc->data; 114 DM dm, cdm; 115 Mat Interp, InterpAdapt; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 /* There is no interpolator for the coarse level */ 120 if (!l) PetscFunctionReturn(0); 121 ierr = KSPGetDM(csmooth, &cdm);CHKERRQ(ierr); 122 ierr = KSPGetDM(fsmooth, &dm);CHKERRQ(ierr); 123 ierr = PCMGGetInterpolation(pc, l, &Interp);CHKERRQ(ierr); 124 125 ierr = DMAdaptInterpolator(cdm, dm, Interp, fsmooth, Nc, fspace, cspace, &InterpAdapt, pc);CHKERRQ(ierr); 126 if (mg->mespMonitor) {ierr = DMCheckInterpolator(dm, InterpAdapt, Nc, cspace, fspace, 0.5/*PETSC_SMALL*/);CHKERRQ(ierr);} 127 ierr = PCMGSetInterpolation(pc, l, InterpAdapt);CHKERRQ(ierr); 128 ierr = PCMGSetRestriction(pc, l, InterpAdapt);CHKERRQ(ierr); 129 ierr = MatDestroy(&InterpAdapt);CHKERRQ(ierr); 130 PetscFunctionReturn(0); 131 } 132 133 /* 134 PCMGRecomputeLevelOperators_Internal - Recomputes Galerkin coarse operator when interpolation is adapted 135 136 Input Parameters: 137 + pc - The PCMG 138 - l - The level l 139 140 Level: developer 141 142 Note: This routine recomputes the Galerkin triple product for the operator on level l. 143 */ 144 PetscErrorCode PCMGRecomputeLevelOperators_Internal(PC pc, PetscInt l) 145 { 146 Mat fA, fB; /* The system and preconditioning operators on level l+1 */ 147 Mat A, B; /* The system and preconditioning operators on level l */ 148 Mat Interp, Restrc; /* The interpolation operator from level l to l+1, and restriction operator from level l+1 to l */ 149 KSP smooth, fsmooth; /* The smoothers on levels l and l+1 */ 150 PCMGGalerkinType galerkin; /* The Galerkin projection flag */ 151 MatReuse reuse = MAT_REUSE_MATRIX; /* The matrices are always assumed to be present already */ 152 PetscBool doA = PETSC_FALSE; /* Updates the system operator */ 153 PetscBool doB = PETSC_FALSE; /* Updates the preconditioning operator (A == B, then update B) */ 154 PetscInt n; /* The number of multigrid levels */ 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = PCMGGetGalerkin(pc, &galerkin);CHKERRQ(ierr); 159 if (galerkin >= PC_MG_GALERKIN_NONE) PetscFunctionReturn(0); 160 ierr = PCMGGetLevels(pc, &n);CHKERRQ(ierr); 161 /* Do not recompute operator for the finest grid */ 162 if (l == n-1) PetscFunctionReturn(0); 163 ierr = PCMGGetSmoother(pc, l, &smooth);CHKERRQ(ierr); 164 ierr = KSPGetOperators(smooth, &A, &B);CHKERRQ(ierr); 165 ierr = PCMGGetSmoother(pc, l+1, &fsmooth);CHKERRQ(ierr); 166 ierr = KSPGetOperators(fsmooth, &fA, &fB);CHKERRQ(ierr); 167 ierr = PCMGGetInterpolation(pc, l+1, &Interp);CHKERRQ(ierr); 168 ierr = PCMGGetRestriction(pc, l+1, &Restrc);CHKERRQ(ierr); 169 if ((galerkin == PC_MG_GALERKIN_PMAT) || (galerkin == PC_MG_GALERKIN_BOTH)) doB = PETSC_TRUE; 170 if ((galerkin == PC_MG_GALERKIN_MAT) || ((galerkin == PC_MG_GALERKIN_BOTH) && (fA != fB))) doA = PETSC_TRUE; 171 if (doA) {ierr = MatGalerkin(Restrc, fA, Interp, reuse, 1.0, &A);CHKERRQ(ierr);} 172 if (doB) {ierr = MatGalerkin(Restrc, fB, Interp, reuse, 1.0, &B);CHKERRQ(ierr);} 173 PetscFunctionReturn(0); 174 } 175