1 static char help[] = 2 "Tests implementation of PetscSpace_Sum by solving the Poisson equations using a PetscSpace_Poly and a PetscSpace_Sum and checking that \ 3 solutions agree up to machine precision.\n\n"; 4 5 #include <petscdmplex.h> 6 #include <petscds.h> 7 #include <petscfe.h> 8 #include <petscsnes.h> 9 /* We are solving the system of equations: 10 * \vec{u} = -\grad{p} 11 * \div{u} = f 12 */ 13 14 /* Exact solutions for linear velocity 15 \vec{u} = \vec{x}; 16 p = -0.5*(\vec{x} \cdot \vec{x}); 17 */ 18 static PetscErrorCode linear_u(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 19 { 20 PetscInt c; 21 22 for (c = 0; c < Nc; ++c) u[c] = x[c]; 23 return 0; 24 } 25 26 static PetscErrorCode linear_p(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 27 { 28 PetscInt d; 29 30 u[0] = 0.; 31 for (d=0; d<dim; ++d) u[0] += -0.5*x[d]*x[d]; 32 return 0; 33 } 34 35 static PetscErrorCode linear_divu(PetscInt dim,PetscReal time,const PetscReal x[],PetscInt Nc,PetscScalar *u,void *ctx) 36 { 37 u[0] = dim; 38 return 0; 39 } 40 41 /* fx_v are the residual functions for the equation \vec{u} = \grad{p}. f0_v is the term <v,u>.*/ 42 static void f0_v(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 43 { 44 PetscInt i; 45 46 for (i=0; i<dim; ++i) f0[i] = u[uOff[0] + i]; 47 } 48 49 /* f1_v is the term <v,-\grad{p}> but we integrate by parts to get <\grad{v}, -p*I> */ 50 static void f1_v(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f1[]) 51 { 52 PetscInt c; 53 54 for (c=0; c<dim; ++c) { 55 PetscInt d; 56 57 for (d=0; d<dim; ++d) f1[c*dim + d] = (c==d) ? -u[uOff[1]] : 0; 58 } 59 } 60 61 /* Residual function for enforcing \div{u} = f. */ 62 static void f0_q_linear(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 63 { 64 PetscScalar rhs,divu=0; 65 PetscInt i; 66 67 (void)linear_divu(dim,t,x,dim,&rhs,NULL); 68 for (i=0; i< dim; ++i) divu += u_x[uOff_x[0]+i*dim+i]; 69 f0[0] = divu-rhs; 70 } 71 72 /* Boundary residual. Dirichlet boundary for u means u_bdy=p*n */ 73 static void f0_bd_u_linear(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],const PetscReal n[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[]) 74 { 75 PetscScalar pressure; 76 PetscInt d; 77 78 (void)linear_p(dim,t,x,dim,&pressure,NULL); 79 for (d=0; d<dim; ++d) f0[d] = pressure*n[d]; 80 } 81 82 /* gx_yz are the jacobian functions obtained by taking the derivative of the y residual w.r.t z*/ 83 static void g0_vu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g0[]) 84 { 85 PetscInt c; 86 87 for (c=0; c<dim; ++c) g0[c*dim + c] = 1.0; 88 } 89 90 static void g1_qu(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g1[]) 91 { 92 PetscInt c; 93 94 for (c=0; c<dim; ++c) g1[c*dim + c] = 1.0; 95 } 96 97 static void g2_vp(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,PetscReal u_tShift,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar g2[]) 98 { 99 PetscInt c; 100 101 for (c=0; c<dim; ++c) g2[c*dim + c] = -1.0; 102 } 103 104 typedef struct 105 { 106 PetscInt dummy; 107 } UserCtx; 108 109 static PetscErrorCode CreateMesh(MPI_Comm comm,UserCtx *user,DM *mesh) 110 { 111 PetscFunctionBegin; 112 PetscCall(DMCreate(comm, mesh)); 113 PetscCall(DMSetType(*mesh, DMPLEX)); 114 PetscCall(DMSetFromOptions(*mesh)); 115 PetscCall(DMSetApplicationContext(*mesh,user)); 116 PetscCall(DMViewFromOptions(*mesh,NULL,"-dm_view")); 117 PetscFunctionReturn(0); 118 } 119 120 /* Setup the system of equations that we wish to solve */ 121 static PetscErrorCode SetupProblem(DM dm,UserCtx *user) 122 { 123 PetscDS ds; 124 DMLabel label; 125 PetscWeakForm wf; 126 const PetscInt id = 1; 127 PetscInt bd; 128 129 PetscFunctionBegin; 130 PetscCall(DMGetDS(dm, &ds)); 131 /* All of these are independent of the user's choice of solution */ 132 PetscCall(PetscDSSetResidual(ds,0,f0_v,f1_v)); 133 PetscCall(PetscDSSetResidual(ds,1,f0_q_linear,NULL)); 134 PetscCall(PetscDSSetJacobian(ds,0,0,g0_vu,NULL,NULL,NULL)); 135 PetscCall(PetscDSSetJacobian(ds,0,1,NULL,NULL,g2_vp,NULL)); 136 PetscCall(PetscDSSetJacobian(ds,1,0,NULL,g1_qu,NULL,NULL)); 137 138 PetscCall(DMGetLabel(dm, "marker", &label)); 139 PetscCall(PetscDSAddBoundary(ds,DM_BC_NATURAL,"Boundary Integral",label,1,&id,0,0,NULL,(void (*)(void))NULL,NULL,user,&bd)); 140 PetscCall(PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 141 PetscCall(PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_linear, 0, NULL)); 142 143 PetscCall(PetscDSSetExactSolution(ds,0,linear_u,NULL)); 144 PetscCall(PetscDSSetExactSolution(ds,1,linear_p,NULL)); 145 PetscFunctionReturn(0); 146 } 147 148 /* Create the finite element spaces we will use for this system */ 149 static PetscErrorCode SetupDiscretization(DM mesh,DM mesh_sum,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 150 { 151 DM cdm = mesh,cdm_sum = mesh_sum; 152 PetscFE u,divu,u_sum,divu_sum; 153 PetscInt dim; 154 PetscBool simplex; 155 156 PetscFunctionBegin; 157 PetscCall(DMGetDimension(mesh, &dim)); 158 PetscCall(DMPlexIsSimplex(mesh, &simplex)); 159 /* Create FE objects and give them names so that options can be set from 160 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */ 161 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,dim,simplex,"velocity_",-1,&u)); 162 PetscCall(PetscObjectSetName((PetscObject)u,"velocity")); 163 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,dim,simplex,"velocity_sum_",-1,&u_sum)); 164 PetscCall(PetscObjectSetName((PetscObject)u_sum,"velocity_sum")); 165 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,1,simplex,"divu_",-1,&divu)); 166 PetscCall(PetscObjectSetName((PetscObject)divu,"divu")); 167 PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,1,simplex,"divu_sum_",-1,&divu_sum)); 168 PetscCall(PetscObjectSetName((PetscObject)divu_sum,"divu_sum")); 169 170 PetscCall(PetscFECopyQuadrature(u,divu)); 171 PetscCall(PetscFECopyQuadrature(u_sum,divu_sum)); 172 173 /* Associate the FE objects with the mesh and setup the system */ 174 PetscCall(DMSetField(mesh,0,NULL,(PetscObject)u)); 175 PetscCall(DMSetField(mesh,1,NULL,(PetscObject)divu)); 176 PetscCall(DMCreateDS(mesh)); 177 PetscCall((*setup)(mesh,user)); 178 179 PetscCall(DMSetField(mesh_sum,0,NULL,(PetscObject)u_sum)); 180 PetscCall(DMSetField(mesh_sum,1,NULL,(PetscObject)divu_sum)); 181 PetscCall(DMCreateDS(mesh_sum)); 182 PetscCall((*setup)(mesh_sum,user)); 183 184 while (cdm) { 185 PetscCall(DMCopyDisc(mesh,cdm)); 186 PetscCall(DMGetCoarseDM(cdm,&cdm)); 187 } 188 189 while (cdm_sum) { 190 PetscCall(DMCopyDisc(mesh_sum,cdm_sum)); 191 PetscCall(DMGetCoarseDM(cdm_sum,&cdm_sum)); 192 } 193 194 /* The Mesh now owns the fields, so we can destroy the FEs created in this 195 * function */ 196 PetscCall(PetscFEDestroy(&u)); 197 PetscCall(PetscFEDestroy(&divu)); 198 PetscCall(PetscFEDestroy(&u_sum)); 199 PetscCall(PetscFEDestroy(&divu_sum)); 200 PetscCall(DMDestroy(&cdm)); 201 PetscCall(DMDestroy(&cdm_sum)); 202 PetscFunctionReturn(0); 203 } 204 205 int main(int argc,char **argv) 206 { 207 UserCtx user; 208 DM dm,dm_sum; 209 SNES snes,snes_sum; 210 Vec u,u_sum; 211 PetscReal errNorm; 212 const PetscReal errTol = PETSC_SMALL; 213 214 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 215 216 /* Set up a snes for the standard approach, one space with 2 components */ 217 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 218 PetscCall(CreateMesh(PETSC_COMM_WORLD,&user,&dm)); 219 PetscCall(SNESSetDM(snes,dm)); 220 221 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */ 222 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes_sum)); 223 PetscCall(CreateMesh(PETSC_COMM_WORLD,&user,&dm_sum)); 224 PetscCall(SNESSetDM(snes_sum,dm_sum)); 225 PetscCall(SetupDiscretization(dm,dm_sum,SetupProblem,&user)); 226 227 /* Set up and solve the system using standard approach. */ 228 PetscCall(DMCreateGlobalVector(dm,&u)); 229 PetscCall(VecSet(u,0.0)); 230 PetscCall(PetscObjectSetName((PetscObject)u,"solution")); 231 PetscCall(DMPlexSetSNESLocalFEM(dm,&user,&user,&user)); 232 PetscCall(SNESSetFromOptions(snes)); 233 PetscCall(DMSNESCheckFromOptions(snes,u)); 234 PetscCall(SNESSolve(snes,NULL,u)); 235 PetscCall(SNESGetSolution(snes,&u)); 236 PetscCall(VecViewFromOptions(u,NULL,"-solution_view")); 237 238 /* Set up and solve the sum space system */ 239 PetscCall(DMCreateGlobalVector(dm_sum,&u_sum)); 240 PetscCall(VecSet(u_sum,0.0)); 241 PetscCall(PetscObjectSetName((PetscObject)u_sum,"solution_sum")); 242 PetscCall(DMPlexSetSNESLocalFEM(dm_sum,&user,&user,&user)); 243 PetscCall(SNESSetFromOptions(snes_sum)); 244 PetscCall(DMSNESCheckFromOptions(snes_sum,u_sum)); 245 PetscCall(SNESSolve(snes_sum,NULL,u_sum)); 246 PetscCall(SNESGetSolution(snes_sum,&u_sum)); 247 PetscCall(VecViewFromOptions(u_sum,NULL,"-solution_sum_view")); 248 249 /* Check if standard solution and sum space solution match to machine precision */ 250 PetscCall(VecAXPY(u_sum,-1,u)); 251 PetscCall(VecNorm(u_sum,NORM_2,&errNorm)); 252 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Sum space provides the same solution as a regular space: %s",(errNorm < errTol) ? "true" : "false")); 253 254 /* Cleanup */ 255 PetscCall(VecDestroy(&u_sum)); 256 PetscCall(VecDestroy(&u)); 257 PetscCall(SNESDestroy(&snes_sum)); 258 PetscCall(SNESDestroy(&snes)); 259 PetscCall(DMDestroy(&dm_sum)); 260 PetscCall(DMDestroy(&dm)); 261 PetscCall(PetscFinalize()); 262 return 0; 263 } 264 265 /*TEST 266 test: 267 suffix: 2d_lagrange 268 requires: triangle 269 args: -velocity_petscspace_degree 1 \ 270 -velocity_petscspace_type poly \ 271 -velocity_petscspace_components 2\ 272 -velocity_petscdualspace_type lagrange \ 273 -divu_petscspace_degree 0 \ 274 -divu_petscspace_type poly \ 275 -divu_petscdualspace_lagrange_continuity false \ 276 -velocity_sum_petscfe_default_quadrature_order 1 \ 277 -velocity_sum_petscspace_degree 1 \ 278 -velocity_sum_petscspace_type sum \ 279 -velocity_sum_petscspace_variables 2 \ 280 -velocity_sum_petscspace_components 2 \ 281 -velocity_sum_petscspace_sum_spaces 2 \ 282 -velocity_sum_petscspace_sum_concatenate true \ 283 -velocity_sum_petscdualspace_type lagrange \ 284 -velocity_sum_sumcomp_0_petscspace_type poly \ 285 -velocity_sum_sumcomp_0_petscspace_degree 1 \ 286 -velocity_sum_sumcomp_0_petscspace_variables 2 \ 287 -velocity_sum_sumcomp_0_petscspace_components 1 \ 288 -velocity_sum_sumcomp_1_petscspace_type poly \ 289 -velocity_sum_sumcomp_1_petscspace_degree 1 \ 290 -velocity_sum_sumcomp_1_petscspace_variables 2 \ 291 -velocity_sum_sumcomp_1_petscspace_components 1 \ 292 -divu_sum_petscspace_degree 0 \ 293 -divu_sum_petscspace_type sum \ 294 -divu_sum_petscspace_variables 2 \ 295 -divu_sum_petscspace_components 1 \ 296 -divu_sum_petscspace_sum_spaces 1 \ 297 -divu_sum_petscspace_sum_concatenate true \ 298 -divu_sum_petscdualspace_lagrange_continuity false \ 299 -divu_sum_sumcomp_0_petscspace_type poly \ 300 -divu_sum_sumcomp_0_petscspace_degree 0 \ 301 -divu_sum_sumcomp_0_petscspace_variables 2 \ 302 -divu_sum_sumcomp_0_petscspace_components 1 \ 303 -dm_refine 0 \ 304 -snes_error_if_not_converged \ 305 -ksp_rtol 1e-10 \ 306 -ksp_error_if_not_converged \ 307 -pc_type fieldsplit\ 308 -pc_fieldsplit_type schur\ 309 -divu_sum_petscdualspace_lagrange_continuity false \ 310 -pc_fieldsplit_schur_precondition full 311 TEST*/ 312