xref: /petsc/src/dm/tutorials/ex21.c (revision 2da392cc7c10228af19ad9843ce5155178acb644) 
1 
2 static char help[] = "DMSwarm-PIC demonstator of advecting points within cell DM defined by a DA or PLEX object \n\
3 Options: \n\
4 -ppcell   : Number of times to sub-divide the reference cell when layout the initial particle coordinates \n\
5 -meshtype : 0 ==> DA , 1 ==> PLEX \n\
6 -nt       : Number of timestep to perform \n\
7 -view     : Write out initial condition and time dependent data \n";
8 
9 #include <petsc.h>
10 #include <petscdm.h>
11 #include <petscdmda.h>
12 #include <petscdmswarm.h>
13 
14 PetscErrorCode pic_advect(PetscInt ppcell,PetscInt meshtype)
15 {
16   PetscErrorCode ierr;
17   const PetscInt dim = 2;
18   DM celldm,swarm;
19   PetscInt tk,nt = 200;
20   PetscBool view = PETSC_FALSE;
21   Vec *pfields;
22   PetscReal minradius;
23   PetscReal dt;
24   PetscReal vel[] = { 1.0, 0.16 };
25   const char *fieldnames[] = { "phi" };
26   PetscViewer viewer;
27 
28   PetscFunctionBegin;
29   ierr = PetscOptionsGetBool(NULL,NULL,"-view",&view,NULL);CHKERRQ(ierr);
30   ierr = PetscOptionsGetInt(NULL,NULL,"-nt",&nt,NULL);CHKERRQ(ierr);
31 
32   /* Create the background cell DM */
33   if (meshtype == 0) { /* DA */
34     PetscInt nxy;
35     PetscInt dof = 1;
36     PetscInt stencil_width = 1;
37 
38     ierr = PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMDA\n");CHKERRQ(ierr);
39     nxy = 33;
40     ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,nxy,nxy,PETSC_DECIDE,PETSC_DECIDE,dof,stencil_width,NULL,NULL,&celldm);CHKERRQ(ierr);
41 
42     ierr = DMDASetElementType(celldm,DMDA_ELEMENT_Q1);CHKERRQ(ierr);
43 
44     ierr = DMSetFromOptions(celldm);CHKERRQ(ierr);
45 
46     ierr = DMSetUp(celldm);CHKERRQ(ierr);
47 
48     ierr = DMDASetUniformCoordinates(celldm,0.0,1.0,0.0,1.0,0.0,1.5);CHKERRQ(ierr);
49 
50     minradius = 1.0/((PetscReal)(nxy-1));
51 
52     ierr = PetscPrintf(PETSC_COMM_WORLD,"DA(minradius) %1.4e\n",(double)minradius);CHKERRQ(ierr);
53   }
54 
55   if (meshtype == 1){ /* PLEX */
56     DM distributedMesh = NULL;
57     PetscInt numComp[] = {1};
58     PetscInt numDof[] = {1,0,0}; /* vert, edge, cell */
59     PetscInt faces[]  = {1,1,1};
60     PetscInt numBC = 0;
61     PetscSection section;
62     Vec cellgeom = NULL;
63     Vec facegeom = NULL;
64 
65     ierr = PetscPrintf(PETSC_COMM_WORLD,"Mesh type: DMPLEX\n");CHKERRQ(ierr);
66     ierr = DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, PETSC_TRUE, faces, NULL, NULL, PETSC_TRUE, &celldm);CHKERRQ(ierr);
67 
68     /* Distribute mesh over processes */
69     ierr = DMPlexDistribute(celldm,0,NULL,&distributedMesh);CHKERRQ(ierr);
70     if (distributedMesh) {
71       ierr = DMDestroy(&celldm);CHKERRQ(ierr);
72       celldm = distributedMesh;
73     }
74 
75     ierr = DMSetFromOptions(celldm);CHKERRQ(ierr);
76 
77     ierr = DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,&section);CHKERRQ(ierr);
78     ierr = DMSetLocalSection(celldm,section);CHKERRQ(ierr);
79 
80     ierr = DMSetUp(celldm);CHKERRQ(ierr);
81 
82     /* Calling DMPlexComputeGeometryFVM() generates the value returned by DMPlexGetMinRadius() */
83     ierr = DMPlexComputeGeometryFVM(celldm,&cellgeom,&facegeom);CHKERRQ(ierr);
84     ierr = DMPlexGetMinRadius(celldm,&minradius);CHKERRQ(ierr);
85     ierr = PetscPrintf(PETSC_COMM_WORLD,"PLEX(minradius) %1.4e\n",(double)minradius);CHKERRQ(ierr);
86     ierr = VecDestroy(&cellgeom);CHKERRQ(ierr);
87     ierr = VecDestroy(&facegeom);CHKERRQ(ierr);
88     ierr = PetscSectionDestroy(&section);CHKERRQ(ierr);
89   }
90 
91   /* Create the DMSwarm */
92   ierr = DMCreate(PETSC_COMM_WORLD,&swarm);CHKERRQ(ierr);
93   ierr = DMSetType(swarm,DMSWARM);CHKERRQ(ierr);
94   ierr = DMSetDimension(swarm,dim);CHKERRQ(ierr);
95 
96   /* Configure swarm to be of type PIC */
97   ierr = DMSwarmSetType(swarm,DMSWARM_PIC);CHKERRQ(ierr);
98   ierr = DMSwarmSetCellDM(swarm,celldm);CHKERRQ(ierr);
99 
100   /* Register two scalar fields within the DMSwarm */
101   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"phi",1,PETSC_REAL);CHKERRQ(ierr);
102   ierr = DMSwarmRegisterPetscDatatypeField(swarm,"region",1,PETSC_REAL);CHKERRQ(ierr);
103   ierr = DMSwarmFinalizeFieldRegister(swarm);CHKERRQ(ierr);
104 
105   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
106   ierr = DMSwarmSetLocalSizes(swarm,4,0);CHKERRQ(ierr);
107 
108   /* Insert swarm coordinates cell-wise */
109   /*ierr = DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,ppcell);CHKERRQ(ierr);*/
110   ierr = DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,ppcell);CHKERRQ(ierr);
111 
112   /* Define initial conditions for th swarm fields "phi" and "region" */
113   {
114     PetscReal *s_coor,*s_phi,*s_region;
115     PetscInt npoints,p;
116 
117     ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr);
118     ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);CHKERRQ(ierr);
119     ierr = DMSwarmGetField(swarm,"phi",NULL,NULL,(void**)&s_phi);CHKERRQ(ierr);
120     ierr = DMSwarmGetField(swarm,"region",NULL,NULL,(void**)&s_region);CHKERRQ(ierr);
121     for (p=0; p<npoints; p++) {
122       PetscReal pos[2];
123       pos[0] = s_coor[2*p+0];
124       pos[1] = s_coor[2*p+1];
125 
126       s_region[p] = 1.0;
127       s_phi[p] = 1.0 + PetscExpReal(-200.0*((pos[0]-0.5)*(pos[0]-0.5) + (pos[1]-0.5)*(pos[1]-0.5)));
128     }
129     ierr = DMSwarmRestoreField(swarm,"region",NULL,NULL,(void**)&s_region);CHKERRQ(ierr);
130     ierr = DMSwarmRestoreField(swarm,"phi",NULL,NULL,(void**)&s_phi);CHKERRQ(ierr);
131     ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);CHKERRQ(ierr);
132   }
133 
134   /* Project initial value of phi onto the mesh */
135   ierr = DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_FALSE);CHKERRQ(ierr);
136 
137   if (view) {
138     /* View swarm all swarm fields using data type PETSC_REAL */
139     ierr = DMSwarmViewXDMF(swarm,"ic_dms.xmf");CHKERRQ(ierr);
140 
141     /* View projected swarm field "phi" */
142     ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
143     ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
144     ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
145     if (meshtype == 0) { /* DA */
146       ierr = PetscViewerFileSetName(viewer,"ic_dmda.vts");CHKERRQ(ierr);
147       ierr = VecView(pfields[0],viewer);CHKERRQ(ierr);
148     }
149     if (meshtype == 1) { /* PLEX */
150       ierr = PetscViewerFileSetName(viewer,"ic_dmplex.vtk");CHKERRQ(ierr);
151       ierr = DMView(celldm,viewer);CHKERRQ(ierr);
152       ierr = VecView(pfields[0],viewer);CHKERRQ(ierr);
153     }
154     ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
155   }
156 
157   ierr = DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
158   ierr = DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
159 
160   dt = 0.5 * minradius / PetscSqrtReal(vel[0]*vel[0] + vel[1]*vel[1]);
161   for (tk=1; tk<=nt; tk++) {
162     PetscReal *s_coor;
163     PetscInt npoints,p;
164 
165     ierr = PetscPrintf(PETSC_COMM_WORLD,"[step %D]\n",tk);CHKERRQ(ierr);
166     /* advect with analytic prescribed (constant) velocity field */
167     ierr = DMSwarmGetLocalSize(swarm,&npoints);CHKERRQ(ierr);
168     ierr = DMSwarmGetField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);CHKERRQ(ierr);
169     for (p=0; p<npoints; p++) {
170       s_coor[2*p+0] += dt * vel[0];
171       s_coor[2*p+1] += dt * vel[1];
172     }
173     ierr = DMSwarmRestoreField(swarm,DMSwarmPICField_coor,NULL,NULL,(void**)&s_coor);CHKERRQ(ierr);
174 
175     ierr = DMSwarmMigrate(swarm,PETSC_TRUE);CHKERRQ(ierr);
176 
177     /* Ad-hoc cell filling algorithm */
178     /*
179        The injection frequency is chosen for default DA case.
180        They will likely not be appropriate for the general case using an unstructure PLEX mesh.
181     */
182     if (tk%10 == 0) {
183       PetscReal dx = 1.0/32.0;
184       PetscInt npoints_dir_x[] = { 32, 1 };
185       PetscReal min[2],max[2];
186 
187       min[0] = 0.5 * dx;  max[0] = 0.5 * dx + 31.0 * dx;
188       min[1] = 0.5 * dx;  max[1] = 0.5 * dx;
189       ierr = DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_x,ADD_VALUES);CHKERRQ(ierr);
190     }
191     if (tk%2 == 0) {
192       PetscReal dx = 1.0/32.0;
193       PetscInt npoints_dir_y[] = { 2, 31 };
194       PetscReal min[2],max[2];
195 
196       min[0] = 0.05 * dx; max[0] = 0.5 * dx;
197       min[1] = 0.5 * dx;  max[1] = 0.5 * dx + 31.0 * dx;
198       ierr = DMSwarmSetPointsUniformCoordinates(swarm,min,max,npoints_dir_y,ADD_VALUES);CHKERRQ(ierr);
199     }
200 
201     /* Project swarm field "phi" onto the cell DM */
202     ierr = DMSwarmProjectFields(swarm,1,fieldnames,&pfields,PETSC_TRUE);CHKERRQ(ierr);
203 
204     if (view) {
205       PetscViewer viewer;
206       char fname[PETSC_MAX_PATH_LEN];
207 
208       /* View swarm fields */
209       ierr = PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dms.xmf",tk);CHKERRQ(ierr);
210       ierr = DMSwarmViewXDMF(swarm,fname);CHKERRQ(ierr);
211 
212       /* View projected field */
213       ierr = PetscViewerCreate(PETSC_COMM_WORLD,&viewer);CHKERRQ(ierr);
214       ierr = PetscViewerSetType(viewer,PETSCVIEWERVTK);CHKERRQ(ierr);
215       ierr = PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
216 
217       if (meshtype == 0) { /* DA */
218         ierr = PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmda.vts",tk);CHKERRQ(ierr);
219         ierr = PetscViewerFileSetName(viewer,fname);CHKERRQ(ierr);
220         ierr = VecView(pfields[0],viewer);CHKERRQ(ierr);
221       }
222       if (meshtype == 1) { /* PLEX */
223         ierr = PetscSNPrintf(fname,PETSC_MAX_PATH_LEN-1,"step%.4D_dmplex.vtk",tk);CHKERRQ(ierr);
224         ierr = PetscViewerFileSetName(viewer,fname);CHKERRQ(ierr);
225         ierr = DMView(celldm,viewer);CHKERRQ(ierr);
226         ierr = VecView(pfields[0],viewer);CHKERRQ(ierr);
227       }
228       ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
229     }
230 
231   }
232   ierr = VecDestroy(&pfields[0]);CHKERRQ(ierr);
233   ierr = PetscFree(pfields);CHKERRQ(ierr);
234   ierr = DMDestroy(&celldm);CHKERRQ(ierr);
235   ierr = DMDestroy(&swarm);CHKERRQ(ierr);
236 
237   PetscFunctionReturn(0);
238 }
239 
240 int main(int argc,char **args)
241 {
242   PetscErrorCode ierr;
243   PetscInt ppcell = 1;
244   PetscInt meshtype = 0;
245 
246   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
247   ierr = PetscOptionsGetInt(NULL,NULL,"-ppcell",&ppcell,NULL);CHKERRQ(ierr);
248   ierr = PetscOptionsGetInt(NULL,NULL,"-meshtype",&meshtype,NULL);CHKERRQ(ierr);
249   if (meshtype > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"-meshtype <value> must be 0 or 1");
250 
251   ierr = pic_advect(ppcell,meshtype);CHKERRQ(ierr);
252 
253   ierr = PetscFinalize();
254   return ierr;
255 }
256