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 PetscBool simplex; 107 PetscInt dim; 108 } UserCtx; 109 110 /* Process command line options and initialize the UserCtx struct */ 111 static PetscErrorCode ProcessOptions(MPI_Comm comm,UserCtx *user) 112 { 113 PetscErrorCode ierr; 114 115 PetscFunctionBegin; 116 /* Default to 2D, triangle mesh.*/ 117 user->simplex = PETSC_TRUE; 118 user->dim = 2; 119 120 ierr = PetscOptionsBegin(comm,"","PetscSpaceSum Tests","PetscSpace");CHKERRQ(ierr); 121 ierr = PetscOptionsBool("-simplex","Whether to use simplices (true) or tensor-product (false) cells in " "the mesh","ex8.c",user->simplex, 122 &user->simplex,NULL);CHKERRQ(ierr); 123 ierr = PetscOptionsInt("-dim","Number of solution dimensions","ex8.c",user->dim,&user->dim,NULL);CHKERRQ(ierr); 124 ierr = PetscOptionsEnd();CHKERRQ(ierr); 125 PetscFunctionReturn(0); 126 } 127 128 static PetscErrorCode CreateMesh(MPI_Comm comm,UserCtx *user,DM *mesh) 129 { 130 PetscErrorCode ierr; 131 DMLabel label; 132 const char *name = "marker"; 133 DM dmDist = NULL; 134 PetscPartitioner part; 135 136 PetscFunctionBegin; 137 /* Create box mesh from user parameters */ 138 ierr = DMPlexCreateBoxMesh(comm,user->dim,user->simplex,NULL,NULL,NULL,NULL,PETSC_TRUE,mesh);CHKERRQ(ierr); 139 140 /* Make sure the mesh gets properly distributed if running in parallel */ 141 ierr = DMPlexGetPartitioner(*mesh,&part);CHKERRQ(ierr); 142 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr); 143 ierr = DMPlexDistribute(*mesh,0,NULL,&dmDist);CHKERRQ(ierr); 144 if (dmDist) { 145 ierr = DMDestroy(mesh);CHKERRQ(ierr); 146 *mesh = dmDist; 147 } 148 149 /* Mark the boundaries, we will need this later when setting up the system of 150 * equations */ 151 ierr = DMCreateLabel(*mesh,name);CHKERRQ(ierr); 152 ierr = DMGetLabel(*mesh,name,&label);CHKERRQ(ierr); 153 ierr = DMPlexMarkBoundaryFaces(*mesh,1,label);CHKERRQ(ierr); 154 ierr = DMPlexLabelComplete(*mesh,label);CHKERRQ(ierr); 155 ierr = DMLocalizeCoordinates(*mesh);CHKERRQ(ierr); 156 ierr = PetscObjectSetName((PetscObject) *mesh,"Mesh");CHKERRQ(ierr); 157 158 /* Get any other mesh options from the command line */ 159 ierr = DMSetApplicationContext(*mesh,user);CHKERRQ(ierr); 160 ierr = DMSetFromOptions(*mesh);CHKERRQ(ierr); 161 ierr = DMViewFromOptions(*mesh,NULL,"-dm_view");CHKERRQ(ierr); 162 163 ierr = DMDestroy(&dmDist);CHKERRQ(ierr); 164 PetscFunctionReturn(0); 165 } 166 167 /* Setup the system of equations that we wish to solve */ 168 static PetscErrorCode SetupProblem(DM dm,UserCtx *user) 169 { 170 PetscDS ds; 171 DMLabel label; 172 PetscWeakForm wf; 173 const PetscInt id = 1; 174 PetscInt bd; 175 PetscErrorCode ierr; 176 177 PetscFunctionBegin; 178 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 179 /* All of these are independent of the user's choice of solution */ 180 ierr = PetscDSSetResidual(ds,0,f0_v,f1_v);CHKERRQ(ierr); 181 ierr = PetscDSSetResidual(ds,1,f0_q_linear,NULL);CHKERRQ(ierr); 182 ierr = PetscDSSetJacobian(ds,0,0,g0_vu,NULL,NULL,NULL);CHKERRQ(ierr); 183 ierr = PetscDSSetJacobian(ds,0,1,NULL,NULL,g2_vp,NULL);CHKERRQ(ierr); 184 ierr = PetscDSSetJacobian(ds,1,0,NULL,g1_qu,NULL,NULL);CHKERRQ(ierr); 185 186 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 187 ierr = PetscDSAddBoundary(ds,DM_BC_NATURAL,"Boundary Integral",label,1,&id,0,0,NULL,(void (*)(void))NULL,NULL,user,&bd);CHKERRQ(ierr); 188 ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 189 ierr = PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, f0_bd_u_linear, 0, NULL);CHKERRQ(ierr); 190 191 ierr = PetscDSSetExactSolution(ds,0,linear_u,NULL);CHKERRQ(ierr); 192 ierr = PetscDSSetExactSolution(ds,1,linear_divu,NULL);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 /* Create the finite element spaces we will use for this system */ 197 static PetscErrorCode SetupDiscretization(DM mesh,DM mesh_sum,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 198 { 199 DM cdm = mesh,cdm_sum = mesh_sum; 200 PetscFE u,divu,u_sum,divu_sum; 201 const PetscInt dim = user->dim; 202 PetscErrorCode ierr; 203 204 PetscFunctionBegin; 205 /* Create FE objects and give them names so that options can be set from 206 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */ 207 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,dim,user->simplex,"velocity_",-1,&u);CHKERRQ(ierr); 208 ierr = PetscObjectSetName((PetscObject)u,"velocity");CHKERRQ(ierr); 209 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,dim,user->simplex,"velocity_sum_",-1,&u_sum);CHKERRQ(ierr); 210 ierr = PetscObjectSetName((PetscObject)u_sum,"velocity_sum");CHKERRQ(ierr); 211 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,1,user->simplex,"divu_",-1,&divu);CHKERRQ(ierr); 212 ierr = PetscObjectSetName((PetscObject)divu,"divu");CHKERRQ(ierr); 213 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,1,user->simplex,"divu_sum_",-1,&divu_sum);CHKERRQ(ierr); 214 ierr = PetscObjectSetName((PetscObject)divu_sum,"divu_sum");CHKERRQ(ierr); 215 216 ierr = PetscFECopyQuadrature(u,divu);CHKERRQ(ierr); 217 ierr = PetscFECopyQuadrature(u_sum,divu_sum);CHKERRQ(ierr); 218 219 /* Associate the FE objects with the mesh and setup the system */ 220 ierr = DMSetField(mesh,0,NULL,(PetscObject)u);CHKERRQ(ierr); 221 ierr = DMSetField(mesh,1,NULL,(PetscObject)divu);CHKERRQ(ierr); 222 ierr = DMCreateDS(mesh);CHKERRQ(ierr); 223 ierr = (*setup)(mesh,user);CHKERRQ(ierr); 224 225 ierr = DMSetField(mesh_sum,0,NULL,(PetscObject)u_sum);CHKERRQ(ierr); 226 ierr = DMSetField(mesh_sum,1,NULL,(PetscObject)divu_sum);CHKERRQ(ierr); 227 ierr = DMCreateDS(mesh_sum);CHKERRQ(ierr); 228 ierr = (*setup)(mesh_sum,user);CHKERRQ(ierr); 229 230 while (cdm) { 231 ierr = DMCopyDisc(mesh,cdm);CHKERRQ(ierr); 232 ierr = DMGetCoarseDM(cdm,&cdm);CHKERRQ(ierr); 233 } 234 235 while (cdm_sum) { 236 ierr = DMCopyDisc(mesh_sum,cdm_sum);CHKERRQ(ierr); 237 ierr = DMGetCoarseDM(cdm_sum,&cdm_sum);CHKERRQ(ierr); 238 } 239 240 /* The Mesh now owns the fields, so we can destroy the FEs created in this 241 * function */ 242 ierr = PetscFEDestroy(&u);CHKERRQ(ierr); 243 ierr = PetscFEDestroy(&divu);CHKERRQ(ierr); 244 ierr = PetscFEDestroy(&u_sum);CHKERRQ(ierr); 245 ierr = PetscFEDestroy(&divu_sum);CHKERRQ(ierr); 246 ierr = DMDestroy(&cdm);CHKERRQ(ierr); 247 ierr = DMDestroy(&cdm_sum);CHKERRQ(ierr); 248 PetscFunctionReturn(0); 249 } 250 251 int main(int argc,char **argv) 252 { 253 UserCtx user; 254 DM dm,dm_sum; 255 SNES snes,snes_sum; 256 Vec u,u_sum; 257 PetscReal errNorm; 258 const PetscReal errTol = PETSC_SMALL; 259 PetscErrorCode ierr; 260 261 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 262 ierr = ProcessOptions(PETSC_COMM_WORLD,&user);CHKERRQ(ierr); 263 264 /* Set up a snes for the standard approach, one space with 2 components */ 265 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 266 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm);CHKERRQ(ierr); 267 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 268 269 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */ 270 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_sum);CHKERRQ(ierr); 271 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm_sum);CHKERRQ(ierr); 272 ierr = SNESSetDM(snes_sum,dm_sum);CHKERRQ(ierr); 273 ierr = SetupDiscretization(dm,dm_sum,SetupProblem,&user);CHKERRQ(ierr); 274 275 /* Set up and solve the system using standard approach. */ 276 ierr = DMCreateGlobalVector(dm,&u);CHKERRQ(ierr); 277 ierr = VecSet(u,0.0);CHKERRQ(ierr); 278 ierr = PetscObjectSetName((PetscObject)u,"solution");CHKERRQ(ierr); 279 ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr); 280 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 281 ierr = DMSNESCheckFromOptions(snes,u);CHKERRQ(ierr); 282 ierr = SNESSolve(snes,NULL,u);CHKERRQ(ierr); 283 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 284 ierr = VecViewFromOptions(u,NULL,"-solution_view");CHKERRQ(ierr); 285 286 /* Set up and solve the sum space system */ 287 ierr = DMCreateGlobalVector(dm_sum,&u_sum);CHKERRQ(ierr); 288 ierr = VecSet(u_sum,0.0);CHKERRQ(ierr); 289 ierr = PetscObjectSetName((PetscObject)u_sum,"solution_sum");CHKERRQ(ierr); 290 ierr = DMPlexSetSNESLocalFEM(dm_sum,&user,&user,&user);CHKERRQ(ierr); 291 ierr = SNESSetFromOptions(snes_sum);CHKERRQ(ierr); 292 ierr = DMSNESCheckFromOptions(snes_sum,u_sum);CHKERRQ(ierr); 293 ierr = SNESSolve(snes_sum,NULL,u_sum);CHKERRQ(ierr); 294 ierr = SNESGetSolution(snes_sum,&u_sum);CHKERRQ(ierr); 295 ierr = VecViewFromOptions(u_sum,NULL,"-solution_sum_view");CHKERRQ(ierr); 296 297 /* Check if standard solution and sum space solution match to machine precision */ 298 ierr = VecAXPY(u_sum,-1,u);CHKERRQ(ierr); 299 ierr = VecNorm(u_sum,NORM_2,&errNorm);CHKERRQ(ierr); 300 ierr = PetscPrintf(PETSC_COMM_WORLD,"Sum space provides the same solution as a regular space: %s",(errNorm < errTol) ? "true" : "false");CHKERRQ( 301 ierr); 302 303 /* Cleanup */ 304 ierr = VecDestroy(&u_sum);CHKERRQ(ierr); 305 ierr = VecDestroy(&u);CHKERRQ(ierr); 306 ierr = SNESDestroy(&snes_sum);CHKERRQ(ierr); 307 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 308 ierr = DMDestroy(&dm_sum);CHKERRQ(ierr); 309 ierr = DMDestroy(&dm);CHKERRQ(ierr); 310 ierr = PetscFinalize(); 311 return ierr; 312 } 313 314 /*TEST 315 test: 316 suffix: 2d_lagrange 317 requires: triangle 318 args: -dim 2 \ 319 -simplex true \ 320 -velocity_petscspace_degree 1 \ 321 -velocity_petscspace_type poly \ 322 -velocity_petscspace_components 2\ 323 -velocity_petscdualspace_type lagrange \ 324 -divu_petscspace_degree 0 \ 325 -divu_petscspace_type poly \ 326 -divu_petscdualspace_lagrange_continuity false \ 327 -velocity_sum_petscfe_default_quadrature_order 1 \ 328 -velocity_sum_petscspace_degree 1 \ 329 -velocity_sum_petscspace_type sum \ 330 -velocity_sum_petscspace_variables 2 \ 331 -velocity_sum_petscspace_components 2 \ 332 -velocity_sum_petscspace_sum_spaces 2 \ 333 -velocity_sum_petscspace_sum_concatenate true \ 334 -velocity_sum_petscdualspace_type lagrange \ 335 -velocity_sum_subspace0_petscspace_type poly \ 336 -velocity_sum_subspace0_petscspace_degree 1 \ 337 -velocity_sum_subspace0_petscspace_variables 2 \ 338 -velocity_sum_subspace0_petscspace_components 1 \ 339 -velocity_sum_subspace1_petscspace_type poly \ 340 -velocity_sum_subspace1_petscspace_degree 1 \ 341 -velocity_sum_subspace1_petscspace_variables 2 \ 342 -velocity_sum_subspace1_petscspace_components 1 \ 343 -divu_sum_petscspace_degree 0 \ 344 -divu_sum_petscspace_type sum \ 345 -divu_sum_petscspace_variables 2 \ 346 -divu_sum_petscspace_components 1 \ 347 -divu_sum_petscspace_sum_spaces 1 \ 348 -divu_sum_petscspace_sum_concatenate true \ 349 -divu_sum_petscdualspace_lagrange_continuity false \ 350 -divu_sum_subspace0_petscspace_type poly \ 351 -divu_sum_subspace0_petscspace_degree 0 \ 352 -divu_sum_subspace0_petscspace_variables 2 \ 353 -divu_sum_subspace0_petscspace_components 1 \ 354 -dm_refine 0 \ 355 -snes_error_if_not_converged \ 356 -ksp_rtol 1e-10 \ 357 -ksp_error_if_not_converged \ 358 -pc_type fieldsplit\ 359 -pc_fieldsplit_type schur\ 360 -divu_sum_petscdualspace_lagrange_continuity false \ 361 -pc_fieldsplit_schur_precondition full 362 TEST*/ 363