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 prob; 171 PetscErrorCode ierr; 172 const PetscInt id=1; 173 174 PetscFunctionBegin; 175 ierr = DMGetDS(dm,&prob);CHKERRQ(ierr); 176 /* All of these are independent of the user's choice of solution */ 177 ierr = PetscDSSetResidual(prob,0,f0_v,f1_v);CHKERRQ(ierr); 178 ierr = PetscDSSetResidual(prob,1,f0_q_linear,NULL);CHKERRQ(ierr); 179 ierr = PetscDSSetJacobian(prob,0,0,g0_vu,NULL,NULL,NULL);CHKERRQ(ierr); 180 ierr = PetscDSSetJacobian(prob,0,1,NULL,NULL,g2_vp,NULL);CHKERRQ(ierr); 181 ierr = PetscDSSetJacobian(prob,1,0,NULL,g1_qu,NULL,NULL);CHKERRQ(ierr); 182 183 ierr = PetscDSAddBoundary(prob,DM_BC_NATURAL,"Boundary Integral","marker",0,0,NULL,(void (*)(void))NULL,NULL,1,&id,user);CHKERRQ(ierr); 184 ierr = PetscDSSetBdResidual(prob,0,f0_bd_u_linear,NULL);CHKERRQ(ierr); 185 ierr = PetscDSSetExactSolution(prob,0,linear_u,NULL);CHKERRQ(ierr); 186 ierr = PetscDSSetExactSolution(prob,1,linear_divu,NULL);CHKERRQ(ierr); 187 PetscFunctionReturn(0); 188 } 189 190 /* Create the finite element spaces we will use for this system */ 191 static PetscErrorCode SetupDiscretization(DM mesh,DM mesh_sum,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 192 { 193 DM cdm = mesh,cdm_sum = mesh_sum; 194 PetscFE u,divu,u_sum,divu_sum; 195 const PetscInt dim = user->dim; 196 PetscErrorCode ierr; 197 198 PetscFunctionBegin; 199 /* Create FE objects and give them names so that options can be set from 200 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */ 201 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,dim,user->simplex,"velocity_",-1,&u);CHKERRQ(ierr); 202 ierr = PetscObjectSetName((PetscObject)u,"velocity");CHKERRQ(ierr); 203 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,dim,user->simplex,"velocity_sum_",-1,&u_sum);CHKERRQ(ierr); 204 ierr = PetscObjectSetName((PetscObject)u_sum,"velocity_sum");CHKERRQ(ierr); 205 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,1,user->simplex,"divu_",-1,&divu);CHKERRQ(ierr); 206 ierr = PetscObjectSetName((PetscObject)divu,"divu");CHKERRQ(ierr); 207 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,1,user->simplex,"divu_sum_",-1,&divu_sum);CHKERRQ(ierr); 208 ierr = PetscObjectSetName((PetscObject)divu_sum,"divu_sum");CHKERRQ(ierr); 209 210 ierr = PetscFECopyQuadrature(u,divu);CHKERRQ(ierr); 211 ierr = PetscFECopyQuadrature(u_sum,divu_sum);CHKERRQ(ierr); 212 213 /* Associate the FE objects with the mesh and setup the system */ 214 ierr = DMSetField(mesh,0,NULL,(PetscObject)u);CHKERRQ(ierr); 215 ierr = DMSetField(mesh,1,NULL,(PetscObject)divu);CHKERRQ(ierr); 216 ierr = DMCreateDS(mesh);CHKERRQ(ierr); 217 ierr = (*setup)(mesh,user);CHKERRQ(ierr); 218 219 ierr = DMSetField(mesh_sum,0,NULL,(PetscObject)u_sum);CHKERRQ(ierr); 220 ierr = DMSetField(mesh_sum,1,NULL,(PetscObject)divu_sum);CHKERRQ(ierr); 221 ierr = DMCreateDS(mesh_sum);CHKERRQ(ierr); 222 ierr = (*setup)(mesh_sum,user);CHKERRQ(ierr); 223 224 while (cdm) { 225 ierr = DMCopyDisc(mesh,cdm);CHKERRQ(ierr); 226 ierr = DMGetCoarseDM(cdm,&cdm);CHKERRQ(ierr); 227 } 228 229 while (cdm_sum) { 230 ierr = DMCopyDisc(mesh_sum,cdm_sum);CHKERRQ(ierr); 231 ierr = DMGetCoarseDM(cdm_sum,&cdm_sum);CHKERRQ(ierr); 232 } 233 234 /* The Mesh now owns the fields, so we can destroy the FEs created in this 235 * function */ 236 ierr = PetscFEDestroy(&u);CHKERRQ(ierr); 237 ierr = PetscFEDestroy(&divu);CHKERRQ(ierr); 238 ierr = PetscFEDestroy(&u_sum);CHKERRQ(ierr); 239 ierr = PetscFEDestroy(&divu_sum);CHKERRQ(ierr); 240 ierr = DMDestroy(&cdm);CHKERRQ(ierr); 241 ierr = DMDestroy(&cdm_sum);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 int main(int argc,char **argv) 246 { 247 UserCtx user; 248 DM dm,dm_sum; 249 SNES snes,snes_sum; 250 Vec u,u_sum; 251 PetscReal errNorm; 252 const PetscReal errTol = PETSC_SMALL; 253 PetscErrorCode ierr; 254 255 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 256 ierr = ProcessOptions(PETSC_COMM_WORLD,&user);CHKERRQ(ierr); 257 258 /* Set up a snes for the standard approach, one space with 2 components */ 259 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 260 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm);CHKERRQ(ierr); 261 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 262 263 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */ 264 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_sum);CHKERRQ(ierr); 265 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm_sum);CHKERRQ(ierr); 266 ierr = SNESSetDM(snes_sum,dm_sum);CHKERRQ(ierr); 267 ierr = SetupDiscretization(dm,dm_sum,SetupProblem,&user);CHKERRQ(ierr); 268 269 /* Set up and solve the system using standard approach. */ 270 ierr = DMCreateGlobalVector(dm,&u);CHKERRQ(ierr); 271 ierr = VecSet(u,0.0);CHKERRQ(ierr); 272 ierr = PetscObjectSetName((PetscObject)u,"solution");CHKERRQ(ierr); 273 ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr); 274 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 275 ierr = DMSNESCheckFromOptions(snes,u);CHKERRQ(ierr); 276 ierr = SNESSolve(snes,NULL,u);CHKERRQ(ierr); 277 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 278 ierr = VecViewFromOptions(u,NULL,"-solution_view");CHKERRQ(ierr); 279 280 /* Set up and solve the sum space system */ 281 ierr = DMCreateGlobalVector(dm_sum,&u_sum);CHKERRQ(ierr); 282 ierr = VecSet(u_sum,0.0);CHKERRQ(ierr); 283 ierr = PetscObjectSetName((PetscObject)u_sum,"solution_sum");CHKERRQ(ierr); 284 ierr = DMPlexSetSNESLocalFEM(dm_sum,&user,&user,&user);CHKERRQ(ierr); 285 ierr = SNESSetFromOptions(snes_sum);CHKERRQ(ierr); 286 ierr = DMSNESCheckFromOptions(snes_sum,u_sum);CHKERRQ(ierr); 287 ierr = SNESSolve(snes_sum,NULL,u_sum);CHKERRQ(ierr); 288 ierr = SNESGetSolution(snes_sum,&u_sum);CHKERRQ(ierr); 289 ierr = VecViewFromOptions(u_sum,NULL,"-solution_sum_view");CHKERRQ(ierr); 290 291 /* Check if standard solution and sum space solution match to machine precision */ 292 ierr = VecAXPY(u_sum,-1,u);CHKERRQ(ierr); 293 ierr = VecNorm(u_sum,NORM_2,&errNorm);CHKERRQ(ierr); 294 ierr = PetscPrintf(PETSC_COMM_WORLD,"Sum space provides the same solution as a regular space: %s",(errNorm < errTol) ? "true" : "false");CHKERRQ( 295 ierr); 296 297 /* Cleanup */ 298 ierr = VecDestroy(&u_sum);CHKERRQ(ierr); 299 ierr = VecDestroy(&u);CHKERRQ(ierr); 300 ierr = SNESDestroy(&snes_sum);CHKERRQ(ierr); 301 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 302 ierr = DMDestroy(&dm_sum);CHKERRQ(ierr); 303 ierr = DMDestroy(&dm);CHKERRQ(ierr); 304 ierr = PetscFinalize(); 305 return ierr; 306 } 307 308 /*TEST 309 test: 310 suffix: 2d_lagrange 311 requires: triangle 312 args: -dim 2 \ 313 -simplex true \ 314 -velocity_petscspace_degree 1 \ 315 -velocity_petscspace_type poly \ 316 -velocity_petscspace_components 2\ 317 -velocity_petscdualspace_type lagrange \ 318 -divu_petscspace_degree 0 \ 319 -divu_petscspace_type poly \ 320 -divu_petscdualspace_lagrange_continuity false \ 321 -velocity_sum_petscfe_default_quadrature_order 1 \ 322 -velocity_sum_petscspace_degree 1 \ 323 -velocity_sum_petscspace_type sum \ 324 -velocity_sum_petscspace_variables 2 \ 325 -velocity_sum_petscspace_components 2 \ 326 -velocity_sum_petscspace_sum_spaces 2 \ 327 -velocity_sum_petscspace_sum_concatenate true \ 328 -velocity_sum_petscdualspace_type lagrange \ 329 -velocity_sum_subspace0_petscspace_type poly \ 330 -velocity_sum_subspace0_petscspace_degree 1 \ 331 -velocity_sum_subspace0_petscspace_variables 2 \ 332 -velocity_sum_subspace0_petscspace_components 1 \ 333 -velocity_sum_subspace1_petscspace_type poly \ 334 -velocity_sum_subspace1_petscspace_degree 1 \ 335 -velocity_sum_subspace1_petscspace_variables 2 \ 336 -velocity_sum_subspace1_petscspace_components 1 \ 337 -divu_sum_petscspace_degree 0 \ 338 -divu_sum_petscspace_type sum \ 339 -divu_sum_petscspace_variables 2 \ 340 -divu_sum_petscspace_components 1 \ 341 -divu_sum_petscspace_sum_spaces 1 \ 342 -divu_sum_petscspace_sum_concatenate true \ 343 -divu_sum_petscdualspace_lagrange_continuity false \ 344 -divu_sum_subspace0_petscspace_type poly \ 345 -divu_sum_subspace0_petscspace_degree 0 \ 346 -divu_sum_subspace0_petscspace_variables 2 \ 347 -divu_sum_subspace0_petscspace_components 1 \ 348 -dm_refine 0 \ 349 -snes_error_if_not_converged \ 350 -ksp_rtol 1e-10 \ 351 -ksp_error_if_not_converged \ 352 -pc_type fieldsplit\ 353 -pc_fieldsplit_type schur\ 354 -divu_sum_petscdualspace_lagrange_continuity false \ 355 -pc_fieldsplit_schur_precondition full 356 TEST*/ 357