Actual source code: ex65dm.c
petsc-3.7.7 2017-09-25
1: static char help[] = "Tests coarsening with DM.\n";
5: #include <petscsys.h>
6: #include <petscvec.h>
7: #include <petscdmda.h>
12: int main(int argc, char **argv)
13: {
14: #if !defined(PETSC_USE_COMPLEX)
16: Vec x,yp1,yp2,yp3,yp4,ym1,ym2,ym3,ym4;
17: PetscReal *values;
18: PetscViewer viewer_in,viewer_outp1,viewer_outp2,viewer_outp3,viewer_outp4;
19: PetscViewer viewer_outm1,viewer_outm2,viewer_outm3,viewer_outm4;
20: DM daf,dac1,dac2,dac3,dac4,daf1,daf2,daf3,daf4;
21: Vec scaling_p1,scaling_p2,scaling_p3,scaling_p4;
22: Mat interp_p1,interp_p2,interp_p3,interp_p4,interp_m1,interp_m2,interp_m3,interp_m4;
23: #endif
25: PetscInitialize(&argc,&argv, (char*)0, help);
26: #if defined(PETSC_USE_COMPLEX)
27: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not for complex numbers");
28: #else
29: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_BOX,1024,1024,PETSC_DECIDE,PETSC_DECIDE, 1, 1,NULL,NULL,&daf);
30: DMCreateGlobalVector(daf,&x);
31: VecGetArray(x,&values);
33: DMCoarsen(daf,PETSC_COMM_WORLD,&dac1);
34: DMCoarsen(dac1,PETSC_COMM_WORLD,&dac2);
35: DMCoarsen(dac2,PETSC_COMM_WORLD,&dac3);
36: DMCoarsen(dac3,PETSC_COMM_WORLD,&dac4);
37: DMRefine(daf,PETSC_COMM_WORLD,&daf1);
38: DMRefine(daf1,PETSC_COMM_WORLD,&daf2);
39: DMRefine(daf2,PETSC_COMM_WORLD,&daf3);
40: DMRefine(daf3,PETSC_COMM_WORLD,&daf4);
42: DMCreateGlobalVector(dac1,&yp1);
43: DMCreateGlobalVector(dac2,&yp2);
44: DMCreateGlobalVector(dac3,&yp3);
45: DMCreateGlobalVector(dac4,&yp4);
46: DMCreateGlobalVector(daf1,&ym1);
47: DMCreateGlobalVector(daf2,&ym2);
48: DMCreateGlobalVector(daf3,&ym3);
49: DMCreateGlobalVector(daf4,&ym4);
51: DMCreateInterpolation(dac1,daf,&interp_p1,&scaling_p1);
52: DMCreateInterpolation(dac2,dac1,&interp_p2,&scaling_p2);
53: DMCreateInterpolation(dac3,dac2,&interp_p3,&scaling_p3);
54: DMCreateInterpolation(dac4,dac3,&interp_p4,&scaling_p4);
55: DMCreateInterpolation(daf,daf1,&interp_m1,NULL);
56: DMCreateInterpolation(daf1,daf2,&interp_m2,NULL);
57: DMCreateInterpolation(daf2,daf3,&interp_m3,NULL);
58: DMCreateInterpolation(daf3,daf4,&interp_m4,NULL);
60: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer_in);
61: PetscViewerBinaryRead(viewer_in,values,1048576,NULL,PETSC_DOUBLE);
62: MatRestrict(interp_p1,x,yp1);
63: VecPointwiseMult(yp1,yp1,scaling_p1);
64: MatRestrict(interp_p2,yp1,yp2);
65: VecPointwiseMult(yp2,yp2,scaling_p2);
66: MatRestrict(interp_p3,yp2,yp3);
67: VecPointwiseMult(yp3,yp3,scaling_p3);
68: MatRestrict(interp_p4,yp3,yp4);
69: VecPointwiseMult(yp4,yp4,scaling_p4);
70: MatRestrict(interp_m1,x,ym1);
71: MatRestrict(interp_m2,ym1,ym2);
72: MatRestrict(interp_m3,ym2,ym3);
73: MatRestrict(interp_m4,ym3,ym4);
75: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_WRITE,&viewer_outp1);
76: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_WRITE,&viewer_outp2);
77: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_WRITE,&viewer_outp3);
78: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_WRITE,&viewer_outp4);
79: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_WRITE,&viewer_outm1);
80: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_WRITE,&viewer_outm2);
81: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_WRITE,&viewer_outm3);
82: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_WRITE,&viewer_outm4);
84: VecView(yp1,viewer_outp1);
85: VecView(x,viewer_outp1);
86: VecView(yp2,viewer_outp2);
87: VecView(yp3,viewer_outp3);
88: VecView(yp4,viewer_outp4);
89: VecView(ym1,viewer_outm1);
90: VecView(ym2,viewer_outm2);
91: VecView(ym3,viewer_outm3);
92: VecView(ym4,viewer_outm4);
94: PetscViewerDestroy(&viewer_in);
95: PetscViewerDestroy(&viewer_outp1);
96: PetscViewerDestroy(&viewer_outp2);
97: PetscViewerDestroy(&viewer_outp3);
98: PetscViewerDestroy(&viewer_outp4);
100: PetscViewerDestroy(&viewer_outm1);
101: PetscViewerDestroy(&viewer_outm2);
102: PetscViewerDestroy(&viewer_outm3);
103: PetscViewerDestroy(&viewer_outm4);
105: VecDestroy(&x);
106: VecDestroy(&yp1);
107: VecDestroy(&yp2);
108: VecDestroy(&yp3);
109: VecDestroy(&yp4);
110: VecDestroy(&ym1);
111: VecDestroy(&ym2);
112: VecDestroy(&ym3);
113: VecDestroy(&ym4);
114: #endif
115: PetscFinalize();
116: return 0;
117: }