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 PetscErrorCode ierr; 112 113 PetscFunctionBegin; 114 ierr = DMCreate(comm, mesh);CHKERRQ(ierr); 115 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr); 116 ierr = DMSetFromOptions(*mesh);CHKERRQ(ierr); 117 ierr = DMSetApplicationContext(*mesh,user);CHKERRQ(ierr); 118 ierr = DMViewFromOptions(*mesh,NULL,"-dm_view");CHKERRQ(ierr); 119 PetscFunctionReturn(0); 120 } 121 122 /* Setup the system of equations that we wish to solve */ 123 static PetscErrorCode SetupProblem(DM dm,UserCtx *user) 124 { 125 PetscDS ds; 126 DMLabel label; 127 PetscWeakForm wf; 128 const PetscInt id = 1; 129 PetscInt bd; 130 PetscErrorCode ierr; 131 132 PetscFunctionBegin; 133 ierr = DMGetDS(dm, &ds);CHKERRQ(ierr); 134 /* All of these are independent of the user's choice of solution */ 135 ierr = PetscDSSetResidual(ds,0,f0_v,f1_v);CHKERRQ(ierr); 136 ierr = PetscDSSetResidual(ds,1,f0_q_linear,NULL);CHKERRQ(ierr); 137 ierr = PetscDSSetJacobian(ds,0,0,g0_vu,NULL,NULL,NULL);CHKERRQ(ierr); 138 ierr = PetscDSSetJacobian(ds,0,1,NULL,NULL,g2_vp,NULL);CHKERRQ(ierr); 139 ierr = PetscDSSetJacobian(ds,1,0,NULL,g1_qu,NULL,NULL);CHKERRQ(ierr); 140 141 ierr = DMGetLabel(dm, "marker", &label);CHKERRQ(ierr); 142 ierr = PetscDSAddBoundary(ds,DM_BC_NATURAL,"Boundary Integral",label,1,&id,0,0,NULL,(void (*)(void))NULL,NULL,user,&bd);CHKERRQ(ierr); 143 ierr = PetscDSGetBoundary(ds, bd, &wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 144 ierr = PetscWeakFormSetIndexBdResidual(wf, label, 1, 0, 0, 0, f0_bd_u_linear, 0, NULL);CHKERRQ(ierr); 145 146 ierr = PetscDSSetExactSolution(ds,0,linear_u,NULL);CHKERRQ(ierr); 147 ierr = PetscDSSetExactSolution(ds,1,linear_divu,NULL);CHKERRQ(ierr); 148 PetscFunctionReturn(0); 149 } 150 151 /* Create the finite element spaces we will use for this system */ 152 static PetscErrorCode SetupDiscretization(DM mesh,DM mesh_sum,PetscErrorCode (*setup)(DM,UserCtx*),UserCtx *user) 153 { 154 DM cdm = mesh,cdm_sum = mesh_sum; 155 PetscFE u,divu,u_sum,divu_sum; 156 PetscInt dim; 157 PetscBool simplex; 158 PetscErrorCode ierr; 159 160 PetscFunctionBegin; 161 ierr = DMGetDimension(mesh, &dim);CHKERRQ(ierr); 162 ierr = DMPlexIsSimplex(mesh, &simplex);CHKERRQ(ierr); 163 /* Create FE objects and give them names so that options can be set from 164 * command line. Each field gets 2 instances (i.e. velocity and velocity_sum)created twice so that we can compare between approaches. */ 165 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,dim,simplex,"velocity_",-1,&u);CHKERRQ(ierr); 166 ierr = PetscObjectSetName((PetscObject)u,"velocity");CHKERRQ(ierr); 167 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,dim,simplex,"velocity_sum_",-1,&u_sum);CHKERRQ(ierr); 168 ierr = PetscObjectSetName((PetscObject)u_sum,"velocity_sum");CHKERRQ(ierr); 169 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh),dim,1,simplex,"divu_",-1,&divu);CHKERRQ(ierr); 170 ierr = PetscObjectSetName((PetscObject)divu,"divu");CHKERRQ(ierr); 171 ierr = PetscFECreateDefault(PetscObjectComm((PetscObject)mesh_sum),dim,1,simplex,"divu_sum_",-1,&divu_sum);CHKERRQ(ierr); 172 ierr = PetscObjectSetName((PetscObject)divu_sum,"divu_sum");CHKERRQ(ierr); 173 174 ierr = PetscFECopyQuadrature(u,divu);CHKERRQ(ierr); 175 ierr = PetscFECopyQuadrature(u_sum,divu_sum);CHKERRQ(ierr); 176 177 /* Associate the FE objects with the mesh and setup the system */ 178 ierr = DMSetField(mesh,0,NULL,(PetscObject)u);CHKERRQ(ierr); 179 ierr = DMSetField(mesh,1,NULL,(PetscObject)divu);CHKERRQ(ierr); 180 ierr = DMCreateDS(mesh);CHKERRQ(ierr); 181 ierr = (*setup)(mesh,user);CHKERRQ(ierr); 182 183 ierr = DMSetField(mesh_sum,0,NULL,(PetscObject)u_sum);CHKERRQ(ierr); 184 ierr = DMSetField(mesh_sum,1,NULL,(PetscObject)divu_sum);CHKERRQ(ierr); 185 ierr = DMCreateDS(mesh_sum);CHKERRQ(ierr); 186 ierr = (*setup)(mesh_sum,user);CHKERRQ(ierr); 187 188 while (cdm) { 189 ierr = DMCopyDisc(mesh,cdm);CHKERRQ(ierr); 190 ierr = DMGetCoarseDM(cdm,&cdm);CHKERRQ(ierr); 191 } 192 193 while (cdm_sum) { 194 ierr = DMCopyDisc(mesh_sum,cdm_sum);CHKERRQ(ierr); 195 ierr = DMGetCoarseDM(cdm_sum,&cdm_sum);CHKERRQ(ierr); 196 } 197 198 /* The Mesh now owns the fields, so we can destroy the FEs created in this 199 * function */ 200 ierr = PetscFEDestroy(&u);CHKERRQ(ierr); 201 ierr = PetscFEDestroy(&divu);CHKERRQ(ierr); 202 ierr = PetscFEDestroy(&u_sum);CHKERRQ(ierr); 203 ierr = PetscFEDestroy(&divu_sum);CHKERRQ(ierr); 204 ierr = DMDestroy(&cdm);CHKERRQ(ierr); 205 ierr = DMDestroy(&cdm_sum);CHKERRQ(ierr); 206 PetscFunctionReturn(0); 207 } 208 209 int main(int argc,char **argv) 210 { 211 UserCtx user; 212 DM dm,dm_sum; 213 SNES snes,snes_sum; 214 Vec u,u_sum; 215 PetscReal errNorm; 216 const PetscReal errTol = PETSC_SMALL; 217 PetscErrorCode ierr; 218 219 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 220 221 /* Set up a snes for the standard approach, one space with 2 components */ 222 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 223 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm);CHKERRQ(ierr); 224 ierr = SNESSetDM(snes,dm);CHKERRQ(ierr); 225 226 /* Set up a snes for the sum space approach, where each subspace of the sum space represents one component */ 227 ierr = SNESCreate(PETSC_COMM_WORLD,&snes_sum);CHKERRQ(ierr); 228 ierr = CreateMesh(PETSC_COMM_WORLD,&user,&dm_sum);CHKERRQ(ierr); 229 ierr = SNESSetDM(snes_sum,dm_sum);CHKERRQ(ierr); 230 ierr = SetupDiscretization(dm,dm_sum,SetupProblem,&user);CHKERRQ(ierr); 231 232 /* Set up and solve the system using standard approach. */ 233 ierr = DMCreateGlobalVector(dm,&u);CHKERRQ(ierr); 234 ierr = VecSet(u,0.0);CHKERRQ(ierr); 235 ierr = PetscObjectSetName((PetscObject)u,"solution");CHKERRQ(ierr); 236 ierr = DMPlexSetSNESLocalFEM(dm,&user,&user,&user);CHKERRQ(ierr); 237 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 238 ierr = DMSNESCheckFromOptions(snes,u);CHKERRQ(ierr); 239 ierr = SNESSolve(snes,NULL,u);CHKERRQ(ierr); 240 ierr = SNESGetSolution(snes,&u);CHKERRQ(ierr); 241 ierr = VecViewFromOptions(u,NULL,"-solution_view");CHKERRQ(ierr); 242 243 /* Set up and solve the sum space system */ 244 ierr = DMCreateGlobalVector(dm_sum,&u_sum);CHKERRQ(ierr); 245 ierr = VecSet(u_sum,0.0);CHKERRQ(ierr); 246 ierr = PetscObjectSetName((PetscObject)u_sum,"solution_sum");CHKERRQ(ierr); 247 ierr = DMPlexSetSNESLocalFEM(dm_sum,&user,&user,&user);CHKERRQ(ierr); 248 ierr = SNESSetFromOptions(snes_sum);CHKERRQ(ierr); 249 ierr = DMSNESCheckFromOptions(snes_sum,u_sum);CHKERRQ(ierr); 250 ierr = SNESSolve(snes_sum,NULL,u_sum);CHKERRQ(ierr); 251 ierr = SNESGetSolution(snes_sum,&u_sum);CHKERRQ(ierr); 252 ierr = VecViewFromOptions(u_sum,NULL,"-solution_sum_view");CHKERRQ(ierr); 253 254 /* Check if standard solution and sum space solution match to machine precision */ 255 ierr = VecAXPY(u_sum,-1,u);CHKERRQ(ierr); 256 ierr = VecNorm(u_sum,NORM_2,&errNorm);CHKERRQ(ierr); 257 ierr = PetscPrintf(PETSC_COMM_WORLD,"Sum space provides the same solution as a regular space: %s",(errNorm < errTol) ? "true" : "false");CHKERRQ( 258 ierr); 259 260 /* Cleanup */ 261 ierr = VecDestroy(&u_sum);CHKERRQ(ierr); 262 ierr = VecDestroy(&u);CHKERRQ(ierr); 263 ierr = SNESDestroy(&snes_sum);CHKERRQ(ierr); 264 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 265 ierr = DMDestroy(&dm_sum);CHKERRQ(ierr); 266 ierr = DMDestroy(&dm);CHKERRQ(ierr); 267 ierr = PetscFinalize(); 268 return ierr; 269 } 270 271 /*TEST 272 test: 273 suffix: 2d_lagrange 274 requires: triangle 275 args: -velocity_petscspace_degree 1 \ 276 -velocity_petscspace_type poly \ 277 -velocity_petscspace_components 2\ 278 -velocity_petscdualspace_type lagrange \ 279 -divu_petscspace_degree 0 \ 280 -divu_petscspace_type poly \ 281 -divu_petscdualspace_lagrange_continuity false \ 282 -velocity_sum_petscfe_default_quadrature_order 1 \ 283 -velocity_sum_petscspace_degree 1 \ 284 -velocity_sum_petscspace_type sum \ 285 -velocity_sum_petscspace_variables 2 \ 286 -velocity_sum_petscspace_components 2 \ 287 -velocity_sum_petscspace_sum_spaces 2 \ 288 -velocity_sum_petscspace_sum_concatenate true \ 289 -velocity_sum_petscdualspace_type lagrange \ 290 -velocity_sum_sumcomp_0_petscspace_type poly \ 291 -velocity_sum_sumcomp_0_petscspace_degree 1 \ 292 -velocity_sum_sumcomp_0_petscspace_variables 2 \ 293 -velocity_sum_sumcomp_0_petscspace_components 1 \ 294 -velocity_sum_sumcomp_1_petscspace_type poly \ 295 -velocity_sum_sumcomp_1_petscspace_degree 1 \ 296 -velocity_sum_sumcomp_1_petscspace_variables 2 \ 297 -velocity_sum_sumcomp_1_petscspace_components 1 \ 298 -divu_sum_petscspace_degree 0 \ 299 -divu_sum_petscspace_type sum \ 300 -divu_sum_petscspace_variables 2 \ 301 -divu_sum_petscspace_components 1 \ 302 -divu_sum_petscspace_sum_spaces 1 \ 303 -divu_sum_petscspace_sum_concatenate true \ 304 -divu_sum_petscdualspace_lagrange_continuity false \ 305 -divu_sum_sumcomp_0_petscspace_type poly \ 306 -divu_sum_sumcomp_0_petscspace_degree 0 \ 307 -divu_sum_sumcomp_0_petscspace_variables 2 \ 308 -divu_sum_sumcomp_0_petscspace_components 1 \ 309 -dm_refine 0 \ 310 -snes_error_if_not_converged \ 311 -ksp_rtol 1e-10 \ 312 -ksp_error_if_not_converged \ 313 -pc_type fieldsplit\ 314 -pc_fieldsplit_type schur\ 315 -divu_sum_petscdualspace_lagrange_continuity false \ 316 -pc_fieldsplit_schur_precondition full 317 TEST*/ 318