#include "petscsnes.h" #include "petscda.h" #include "petscdmmg.h" /* User-defined routines and data structures */ typedef struct { PetscScalar u,v,phi; } Field; extern PetscErrorCode FormInitialGuess(DMMG,Vec); extern PetscErrorCode FormFunctionLocal(DALocalInfo*,Field**,Field**,void*); extern PetscErrorCode VecView_VTK(Vec,const char []); typedef struct { PassiveReal lidvelocity,prandtl,xi; /* physical parameters */ PetscTruth draw_contours; /* flag - 1 indicates drawing contours */ } AppCtx; #undef __FUNCT__ #define __FUNCT__ "main" int main(int argc,char **argv) { DMMG *dmmg; /* multilevel grid structure */ AppCtx user; /* user-defined work context */ PetscInt mx,my,its,iter,nlevels=5; PetscErrorCode ierr; MPI_Comm comm; SNES snes; DA da; PetscViewer viewer; Vec sol,zerosol; VecScatter ctx; ierr = PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);if (ierr) return(1); comm = PETSC_COMM_WORLD; ierr = PetscOptionsGetInt(PETSC_NULL,"-nlevels",&nlevels,PETSC_NULL);CHKERRQ(ierr); PreLoadBegin(PETSC_TRUE,"SetUp"); ierr = DMMGCreate(comm,nlevels,&user,&dmmg);CHKERRQ(ierr); /* Create distributed array multigrid object (DMMG) to manage parallel grid and vectors for principal unknowns (x) and governing residuals (f) */ ierr = DACreate2d(comm,DA_NONPERIODIC,DA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,3,1,0,0,&da);CHKERRQ(ierr); /* /\* is this the coarsest grid*\/ */ ierr = DMMGSetDM(dmmg,(DM)da);CHKERRQ(ierr); ierr = DADestroy(da);CHKERRQ(ierr); ierr = DAGetInfo(DMMGGetDA(dmmg),0,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); /* /\*0 for dimension of the DA and also ignores DOF. Calling to get mx and my only. *\/ */ /* Problem parameters (velocity of lid, prandtl, and xi numbers) */ user.lidvelocity = 1.0/(mx*my); user.prandtl = 1.0; user.xi = 1.0; ierr = DASetFieldName(DMMGGetDA(dmmg),0,"x-velocity");CHKERRQ(ierr); ierr = DASetFieldName(DMMGGetDA(dmmg),1,"y-velocity");CHKERRQ(ierr); ierr = DASetFieldName(DMMGGetDA(dmmg),2,"Phi");CHKERRQ(ierr); ierr = DMMGSetSNESLocal(dmmg,FormFunctionLocal,0,0,0);CHKERRQ(ierr); ierr = DMMGSetFromOptions(dmmg);CHKERRQ(ierr); ierr = PetscPrintf(comm,"lid velocity = %G, prandtl # = %G, xi # = %G\n", user.lidvelocity,user.prandtl,user.xi);CHKERRQ(ierr); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Solve the nonlinear system - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMMGSetInitialGuess(dmmg,FormInitialGuess);CHKERRQ(ierr); PreLoadStage("Solve"); ierr = DMMGSolve(dmmg);CHKERRQ(ierr); snes = DMMGGetSNES(dmmg); ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); ierr = PetscPrintf(comm,"Number of Newton iterations = %D\n", its);CHKERRQ(ierr); ierr = VecView_VTK(DMMGGetx(dmmg), "solution.vtk");CHKERRQ(ierr); /* ierr= PetscViewerASCIIOpen(PETSC_COMM_WORLD,"solution.vtk",&viewer);CHKERRQ(ierr); */ /* ierr = VecView(DMMGGetx(dmmg), viewer);CHKERRQ(ierr); */ /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Free work space. All PETSc objects should be destroyed when they are no longer needed. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ ierr = DMMGDestroy(dmmg);CHKERRQ(ierr); PreLoadEnd(); ierr = PetscFinalize();CHKERRQ(ierr); return 0; } /* ------------------------------------------------------------------- */ #undef __FUNCT__ #define __FUNCT__ "FormInitialGuess" /* FormInitialGuess - Forms initial approximation. Input Parameters: user - user-defined application context X - vector Output Parameter: X - vector */ PetscErrorCode FormInitialGuess(DMMG dmmg,Vec X) { AppCtx *user = (AppCtx*)dmmg->user; DA da = (DA)dmmg->dm; PetscInt i,j,mx,xs,ys,xm,ym; PetscErrorCode ierr; PetscReal xi,dx; Field **x; xi = user->xi; ierr = DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);CHKERRQ(ierr); dx = 1.0/(mx-1); /* Get local grid boundaries (for 2-dimensional DA): xs, ys - starting grid indices (no ghost points) xm, ym - widths of local grid (no ghost points) */ ierr = DAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);CHKERRQ(ierr); /* Get a pointer to vector data. - For default PETSc vectors, VecGetArray() returns a pointer to the data array. Otherwise, the routine is implementation dependent. - You MUST call VecRestoreArray() when you no longer need access to the array. */ ierr = DAVecGetArray(da,X,&x);CHKERRQ(ierr); /* Compute initial guess over the locally owned part of the grid Initial condition is motionless fluid and equilibrium temperature */ for (j=ys; j0)*i*dx; /* x[j][i].temp = (xi>0)*i*dx; */ } } /* Restore vector */ ierr = DAVecRestoreArray(da,X,&x);CHKERRQ(ierr); return 0; } PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr) /* x is the input and f is the output, they are both double pointers */ { AppCtx *user = (AppCtx*)ptr; /* Some form of pointer to pointer? */ PetscErrorCode ierr; PetscInt xints,xinte,yints,yinte,i,j; PetscReal hx,hy,dhx,dhy,hxdhy,hydhx; PetscReal xi,prandtl,lid; PetscScalar u,uxx,uyy,vx,vy,avx,avy,vxp,vxm,vyp,vym; PetscScalar phir,phil,phin,phis; PetscFunctionBegin; xi = user->xi; prandtl = user->prandtl; lid = user->lidvelocity; /* because user is not a direct appctx we use -> instead of . */ /* Define mesh intervals ratios for uniform grid. Note: FD formulae below are normalized by multiplying through by local volume element (i.e. hx*hy) to obtain coefficients O(1) in two dimensions. */ dhx = (PetscReal)(info->mx-1); dhy = (PetscReal)(info->my-1); /* info is a pointer hence -> instead of . */ hx = 1.0/dhx; hy = 1.0/dhy; hxdhy = hx*dhy; hydhx = hy*dhx; xints = info->xs; xinte = info->xs+info->xm; yints = info->ys; yinte = info->ys+info->ym; /* Test whether we are on the bottom edge of the global array */ if (yints == 0) { j = 0; yints = yints + 1; /* bottom edge */ for (i=info->xs; ixs+info->xm; i++) { f[j][i].u = x[j][i].u; f[j][i].v = x[j][i].v; f[j][i].phi = x[j][i].phi;/* + (x[j+1][i].u - x[j][i].u)*dhy; */ /* f[j][i].u=0.0; */ /* f[j][i].v=0.0; */ /* f[j][i].phi=0.1; */ } } /* Test whether we are on the top edge of the global array */ if (yinte == info->my) { j = info->my - 1; yinte = yinte - 1; /* top edge */ for (i=info->xs; ixs+info->xm; i++) { f[j][i].u = x[j][i].u - lid; f[j][i].v = x[j][i].v; f[j][i].phi = x[j][i].phi;/* + (x[j][i].u - x[j-1][i].u)*dhy; */ /* f[j][i].u=-lid; */ /* f[j][i].v=0.0; */ /* f[j][i].phi=0.1; */ } } /* Test whether we are on the left edge of the global array */ if (xints == 0) { i = 0; xints = xints + 1; /* left edge */ for (j=info->ys; jys+info->ym; j++) { f[j][i].u = x[j][i].u; f[j][i].v = x[j][i].v; f[j][i].phi = x[j][i].phi;/* - (x[j][i+1].v - x[j][i].v)*dhx; */ } } /* Test whether we are on the right edge of the global array */ if (xinte == info->mx) { i = info->mx - 1; xinte = xinte - 1; /* right edge */ for (j=info->ys; jys+info->ym; j++) { f[j][i].u = x[j][i].u; f[j][i].v = x[j][i].v; f[j][i].phi = x[j][i].phi;/* - (x[j][i].v - x[j][i-1].v)*dhx; */ } } /* Compute over the interior points */ for (j=yints; jym*info->xm);CHKERRQ(ierr); PetscFunctionReturn(0); } #undef __FUNCT__ #define __FUNCT__ "VecView_VTK" PetscErrorCode VecView_VTK(Vec x, const char filename[]) { MPI_Comm comm; DA da; Vec coords; PetscViewer viewer; PetscScalar *array, *values; PetscInt n, N, maxn, mx, my, dof; PetscInt i, p; MPI_Status status; PetscMPIInt rank, size, tag, nn; PetscErrorCode ierr; PetscFunctionBegin; ierr = PetscObjectGetComm((PetscObject) x, &comm);CHKERRQ(ierr); ierr = PetscViewerASCIIOpen(comm, filename, &viewer);CHKERRQ(ierr); ierr = VecGetSize(x, &N); CHKERRQ(ierr); ierr = VecGetLocalSize(x, &n); CHKERRQ(ierr); ierr = PetscObjectQuery((PetscObject) x, "DA", (PetscObject *) &da);CHKERRQ(ierr); if (!da) SETERRQ(PETSC_ERR_ARG_WRONG,"Vector not generated from a DA"); ierr = DAGetInfo(da, 0, &mx, &my, 0,0,0,0, &dof,0,0,0);CHKERRQ(ierr); /* getting information on the global x y size*/ ierr = PetscViewerASCIIPrintf(viewer, "# vtk DataFile Version 2.0\n");CHKERRQ(ierr); /* ierr = PetscViewerASCIIPrintf(viewer, "Inhomogeneous Poisson Equation with %s boundary conditions\n", bcName);CHKERRQ(ierr);*/ ierr = PetscViewerASCIIPrintf(viewer, "ASCII\n");CHKERRQ(ierr); /* get coordinates of nodes */ ierr = DAGetCoordinates(da, &coords);CHKERRQ(ierr); if (!coords) { ierr = DASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, PETSC_NULL,PETSC_NULL );CHKERRQ(ierr); ierr = DAGetCoordinates(da, &coords);CHKERRQ(ierr); } /* If not, set uniform nodes x[0,1], y[0,1]*/ ierr = PetscViewerASCIIPrintf(viewer, "DATASET RECTILINEAR_GRID\n");CHKERRQ(ierr); /* ierr = PetscViewerASCIIPrintf(viewer, "DATASET STRUCTURED_POINTS\n");CHKERRQ(ierr); */ ierr = PetscViewerASCIIPrintf(viewer, "DIMENSIONS %d %d %d\n", mx, my, 1);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "ORIGIN %d %d %d\n", 0, 0, 0);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "SPACING %d %d %d\n", 1, 1, 1);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", mx*my);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", dof);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");CHKERRQ(ierr); ierr = VecGetArray(coords, &array);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "X_COORDINATES %d double\n", mx);CHKERRQ(ierr); for(i = 0; i < mx; i++) { ierr = PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));CHKERRQ(ierr); } /* for(i = 0; i < mx*my-1; i++) { */ /* ierr = PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));CHKERRQ(ierr); */ /* } */ ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "Y_COORDINATES %d double\n", my);CHKERRQ(ierr); for(i = 0; i < my; i++) { ierr = PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*mx*2+1]));CHKERRQ(ierr); } /* for(i = 0; i < mx*my-1; i++) { */ /* ierr = PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2+1]));CHKERRQ(ierr); */ /* } */ ierr = PetscViewerASCIIPrintf(viewer, "\n");CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "Z_COORDINATES %d double\n", 1);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "%G\n", 0.0);CHKERRQ(ierr); ierr = VecRestoreArray(coords, &array);CHKERRQ(ierr); ierr = VecDestroy(coords);CHKERRQ(ierr); ierr = VecGetArray(x, &array);CHKERRQ(ierr); /* Determine maximum message to arrive */ ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); ierr = MPI_Reduce(&n, &maxn, 1, MPIU_INT, MPI_MAX, 0, comm);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "rank %d\n", rank);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "size %d\n", size);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "maxn %d\n", maxn);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "Local size of x, n %d\n", n);CHKERRQ(ierr); ierr = PetscViewerASCIIPrintf(viewer, "Global size of x, N %d\n", N);CHKERRQ(ierr); tag = ((PetscObject) viewer)->tag; if (!rank) { ierr = PetscMalloc((maxn+1) * sizeof(PetscScalar), &values);CHKERRQ(ierr); for(i = 0; i < n; i++) { ierr = PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));CHKERRQ(ierr); } for(p = 1; p < size; p++) { ierr = MPI_Recv(values, (PetscMPIInt) n, MPIU_SCALAR, p, tag, comm, &status);CHKERRQ(ierr); ierr = MPI_Get_count(&status, MPIU_SCALAR, &nn);CHKERRQ(ierr); for(i = 0; i < nn; i++) { ierr = PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));CHKERRQ(ierr); } } ierr = PetscFree(values);CHKERRQ(ierr); } else { ierr = MPI_Send(array, n, MPIU_SCALAR, 0, tag, comm);CHKERRQ(ierr); } ierr = VecRestoreArray(x, &array);CHKERRQ(ierr); ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr); PetscFunctionReturn(0); }