Actual source code: ex65dm.c

petsc-master 2016-12-09
Report Typos and Errors
  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: {
 15: #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);if (ierr) return ierr;
 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:   DMSetFromOptions(daf);
 31:   DMSetUp(daf);
 32:   DMCreateGlobalVector(daf,&x);
 33:   VecGetArray(x,&values);

 35:   DMCoarsen(daf,PETSC_COMM_WORLD,&dac1);
 36:   DMCoarsen(dac1,PETSC_COMM_WORLD,&dac2);
 37:   DMCoarsen(dac2,PETSC_COMM_WORLD,&dac3);
 38:   DMCoarsen(dac3,PETSC_COMM_WORLD,&dac4);
 39:   DMRefine(daf,PETSC_COMM_WORLD,&daf1);
 40:   DMRefine(daf1,PETSC_COMM_WORLD,&daf2);
 41:   DMRefine(daf2,PETSC_COMM_WORLD,&daf3);
 42:   DMRefine(daf3,PETSC_COMM_WORLD,&daf4);

 44:   DMCreateGlobalVector(dac1,&yp1);
 45:   DMCreateGlobalVector(dac2,&yp2);
 46:   DMCreateGlobalVector(dac3,&yp3);
 47:   DMCreateGlobalVector(dac4,&yp4);
 48:   DMCreateGlobalVector(daf1,&ym1);
 49:   DMCreateGlobalVector(daf2,&ym2);
 50:   DMCreateGlobalVector(daf3,&ym3);
 51:   DMCreateGlobalVector(daf4,&ym4);

 53:   DMCreateInterpolation(dac1,daf,&interp_p1,&scaling_p1);
 54:   DMCreateInterpolation(dac2,dac1,&interp_p2,&scaling_p2);
 55:   DMCreateInterpolation(dac3,dac2,&interp_p3,&scaling_p3);
 56:   DMCreateInterpolation(dac4,dac3,&interp_p4,&scaling_p4);
 57:   DMCreateInterpolation(daf,daf1,&interp_m1,NULL);
 58:   DMCreateInterpolation(daf1,daf2,&interp_m2,NULL);
 59:   DMCreateInterpolation(daf2,daf3,&interp_m3,NULL);
 60:   DMCreateInterpolation(daf3,daf4,&interp_m4,NULL);

 62:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi",FILE_MODE_READ,&viewer_in);
 63:   PetscViewerBinaryRead(viewer_in,values,1048576,NULL,PETSC_DOUBLE);
 64:   MatRestrict(interp_p1,x,yp1);
 65:   VecPointwiseMult(yp1,yp1,scaling_p1);
 66:   MatRestrict(interp_p2,yp1,yp2);
 67:   VecPointwiseMult(yp2,yp2,scaling_p2);
 68:   MatRestrict(interp_p3,yp2,yp3);
 69:   VecPointwiseMult(yp3,yp3,scaling_p3);
 70:   MatRestrict(interp_p4,yp3,yp4);
 71:   VecPointwiseMult(yp4,yp4,scaling_p4);
 72:   MatRestrict(interp_m1,x,ym1);
 73:   MatRestrict(interp_m2,ym1,ym2);
 74:   MatRestrict(interp_m3,ym2,ym3);
 75:   MatRestrict(interp_m4,ym3,ym4);

 77:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi1",FILE_MODE_WRITE,&viewer_outp1);
 78:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi2",FILE_MODE_WRITE,&viewer_outp2);
 79:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi3",FILE_MODE_WRITE,&viewer_outp3);
 80:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phi4",FILE_MODE_WRITE,&viewer_outp4);
 81:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim1",FILE_MODE_WRITE,&viewer_outm1);
 82:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim2",FILE_MODE_WRITE,&viewer_outm2);
 83:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim3",FILE_MODE_WRITE,&viewer_outm3);
 84:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,"phim4",FILE_MODE_WRITE,&viewer_outm4);

 86:   VecView(yp1,viewer_outp1);
 87:   VecView(x,viewer_outp1);
 88:   VecView(yp2,viewer_outp2);
 89:   VecView(yp3,viewer_outp3);
 90:   VecView(yp4,viewer_outp4);
 91:   VecView(ym1,viewer_outm1);
 92:   VecView(ym2,viewer_outm2);
 93:   VecView(ym3,viewer_outm3);
 94:   VecView(ym4,viewer_outm4);

 96:   PetscViewerDestroy(&viewer_in);
 97:   PetscViewerDestroy(&viewer_outp1);
 98:   PetscViewerDestroy(&viewer_outp2);
 99:   PetscViewerDestroy(&viewer_outp3);
100:   PetscViewerDestroy(&viewer_outp4);

102:   PetscViewerDestroy(&viewer_outm1);
103:   PetscViewerDestroy(&viewer_outm2);
104:   PetscViewerDestroy(&viewer_outm3);
105:   PetscViewerDestroy(&viewer_outm4);

107:   VecDestroy(&x);
108:   VecDestroy(&yp1);
109:   VecDestroy(&yp2);
110:   VecDestroy(&yp3);
111:   VecDestroy(&yp4);
112:   VecDestroy(&ym1);
113:   VecDestroy(&ym2);
114:   VecDestroy(&ym3);
115:   VecDestroy(&ym4);
116: #endif
117:   PetscFinalize();
118:   return ierr;
119: }