[cig-commits] commit: added grid.dx grid.dz to all functions (including IO and Vis)
Mercurial
hg at geodynamics.org
Thu May 29 11:26:22 PDT 2008
changeset: 127:4adadc50962c
user: Marc Spiegelman <mspieg at ldeo.columbia.edu>
date: Sun Mar 16 23:47:55 2008 -0400
files: matlab/viewSolitaryWave2D.m matlab/viewSolitaryWave2DSurf.m src/FieldFunctions.c src/Params.c src/PetscSimOutput.c src/SolitaryWave2DSL.c src/SolitaryWave2DSL.h src/makeOutput.c
description:
added grid.dx grid.dz to all functions (including IO and Vis)
M matlab/viewSolitaryWave2D.m
M matlab/viewSolitaryWave2DSurf.m
M src/FieldFunctions.c
M src/Params.c
M src/SolitaryWave2DSL.c
M src/SolitaryWave2DSL.h
M src/makeOutput.c
R src/PetscSimOutput.c
diff -r b68f2b824135 -r 4adadc50962c matlab/viewSolitaryWave2D.m
--- a/matlab/viewSolitaryWave2D.m Sun Mar 16 22:44:05 2008 -0400
+++ b/matlab/viewSolitaryWave2D.m Sun Mar 16 23:47:55 2008 -0400
@@ -14,8 +14,8 @@ function Set = viewSolitaryWave2D(fileNa
% calculate the x and z coordinates and add to Set
[m,n]=size(Set.field.phi);
- x = linspace(0,Set.par.width,n);
- y = linspace(0,Set.par.height,m);
+ x = Set.par.dx*[0:n-1]';
+ y = Set.par.dz*[0:m-1]';
Set.x = x;
Set.y = y;
diff -r b68f2b824135 -r 4adadc50962c matlab/viewSolitaryWave2DSurf.m
--- a/matlab/viewSolitaryWave2DSurf.m Sun Mar 16 22:44:05 2008 -0400
+++ b/matlab/viewSolitaryWave2DSurf.m Sun Mar 16 23:47:55 2008 -0400
@@ -15,8 +15,8 @@ function Set = viewSolitaryWave2DSurf(fi
% calculate the x and z coordinates and add to Set
[m,n]=size(Set.field.phi);
- x = linspace(0,Set.par.width,n);
- y = linspace(0,Set.par.height,m);
+ x = Set.par.dx*[0:n-1]';
+ y = Set.par.dz*[0:m-1]';
Set.x = x;
Set.y = y;
diff -r b68f2b824135 -r 4adadc50962c src/FieldFunctions.c
--- a/src/FieldFunctions.c Sun Mar 16 22:44:05 2008 -0400
+++ b/src/FieldFunctions.c Sun Mar 16 23:47:55 2008 -0400
@@ -41,8 +41,8 @@ int Initialize(DMMG *dmmg)
/* Grid parameters */
da = grid->daC;
ierr = DAGetInfo(da,0,&ni,&nj,0,0,0,0,0,0,0,0); CHKERRQ(ierr);
- dx = param->width/ni;
- dz = param->height/nj;
+ dx = grid->dx;
+ dz = grid->dz;
if (param->inputFromFile) { /* Load initial field from a file */
ierr = PetscViewerBinaryOpen(comm,param->inputFilename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
@@ -174,12 +174,13 @@ PetscErrorCode FormFunctionLocal(DALocal
PetscErrorCode FormFunctionLocal(DALocalInfo *info,Field **x,Field **f,void *ptr)
/*---------------------------------------------------------------------*/
{
- AppCtx *user = (AppCtx*) ptr;
+ AppCtx *user = (AppCtx*) ptr;
PetscBag bag = user->bag;
Parameter *param;
- TStepCtx *tsCtx = user->tsCtx;
- Field **xstar;
- DA da=info->da;
+ TStepCtx *tsCtx = user->tsCtx;
+ GridInfo *grid = user->grid;
+ Field **xstar;
+ DA da=info->da;
int i,j;
int ierr, mx, my;
@@ -199,18 +200,18 @@ PetscErrorCode FormFunctionLocal(DALocal
mx = info->mx;
my = info->my;
- dx = param->width/(PetscReal) mx;
- dz = param->height/(PetscReal) my;
-
- dzdx2=dz*dz/dx/dx; /* grid aspect ratio squared */
- dz2=dz*dz;
+ dx = grid->dx;
+ dz = grid->dz;
+
+ dzdx2 = dz*dz/dx/dx; /* grid aspect ratio squared */
+ dz2 = dz*dz;
xints = info->xs; xinte = info->xs+info->xm;
yints = info->ys; yinte = info->ys+info->ym;
/* get time stepping parameters */
- t=tsCtx->t;
- dt=tsCtx->dt;
+ t = tsCtx->t;
+ dt = tsCtx->dt;
/* if a full time step set time-step constant */
if ( tsCtx->takeFullStep) {
@@ -259,24 +260,25 @@ PetscErrorCode FormFunctionLocal(DALocal
#define __FUNCT__ "FormJacobianLocal"
PetscErrorCode FormJacobianLocal(DALocalInfo *info,Field **x,Mat jac,void *ptr)
{
- AppCtx *user = (AppCtx*) ptr;
- PetscBag bag = user->bag;
- Parameter *param;
- TStepCtx *tsCtx = user->tsCtx;
+ AppCtx *user = (AppCtx*) ptr;
+ PetscBag bag = user->bag;
+ Parameter *param;
+ TStepCtx *tsCtx = user->tsCtx;
+ GridInfo *grid = user->grid;
MatStencil col[10],row;
PetscScalar v[10];
int nColumns;
- int i,j,ierr;
- int xints,xinte,yints,yinte;
- PetscReal nExp;
- int mx,my;
- PetscScalar phlfE,phlfW,phlfN,phlfS; /* porosity at the half-points */
- PetscScalar khlfE,khlfW,khlfN,khlfS,sum; /* permeability at the half-points */
- PetscScalar dx, dz, dzdx2, dz2;
- PetscScalar t,dt,halfDt,halfNExp;
- PetscScalar one = 1.0, cmpC;
+ int i,j,ierr;
+ int xints,xinte,yints,yinte;
+ PetscReal nExp;
+ int mx,my;
+ PetscScalar phlfE,phlfW,phlfN,phlfS; /* porosity at the half-points */
+ PetscScalar khlfE,khlfW,khlfN,khlfS,sum; /* permeability at the half-points */
+ PetscScalar dx, dz, dzdx2, dz2;
+ PetscScalar t,dt,halfDt,halfNExp;
+ PetscScalar one = 1.0, cmpC;
PetscFunctionBegin;
@@ -289,10 +291,10 @@ PetscErrorCode FormJacobianLocal(DALocal
my = info->my;
/* Define global and local grid parameters */
- dx = param->width/(PetscReal) (mx); /* for doubly periodic domains if period is width then there are mx+1 grid points per width */
- dz = param->height/(PetscReal) (my); /* ditto */
- dzdx2=dz*dz/dx/dx; /* grid aspect ratio squared */
- dz2=dz*dz;
+ dx = grid->dx;
+ dz = grid->dz;
+ dzdx2 = dz*dz/dx/dx; /* grid aspect ratio squared */
+ dz2 = dz*dz;
xints = info->xs; xinte = info->xs+info->xm;
yints = info->ys; yinte = info->ys+info->ym;
@@ -373,125 +375,6 @@ PetscErrorCode FormJacobianLocal(DALocal
MatSetValuesStencil(jac,1,&row,nColumns,col,v,INSERT_VALUES);
- }
- }
-
- MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
- MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
- PetscFunctionReturn(0);
-}
-
-
-/*
- ---------------------------------------------------------------------
- Calculates the local part of the approximate Jacobian neglecting the B block
- arising from Newton on dF/d\phi
- ---------------------------------------------------------------------
-*/
-#undef __FUNCT__
-#define __FUNCT__ "FormJacobianLocalApprox"
-PetscErrorCode FormJacobianLocalApprox(DALocalInfo *info,Field **x,Mat jac,void *ptr)
-{
- AppCtx *user = (AppCtx*) ptr;
- PetscBag bag = user->bag;
- Parameter *param;
- TStepCtx *tsCtx = user->tsCtx;
- MatStencil col[5],row;
- PetscScalar v[5];
-
-
- int i,j;
- int ierr;
- int xints,xinte,yints,yinte;
- PetscReal nExp;
- int mx,my;
- PetscScalar phlfE,phlfW,phlfN,phlfS; /* porosity at the half-points */
- PetscScalar khlfE,khlfW,khlfN,khlfS,sum; /* permeability at the half-points */
- PetscScalar dx, dz, dzdx2, dz2;
- PetscScalar t,dt,halfDt;
- PetscScalar one = 1.0;
-
- PetscFunctionBegin;
-
- /* retrieve parameters from the bag */
- ierr = PetscBagGetData(bag,(void**)¶m);CHKERRQ(ierr);
- nExp = param->nExp;
-
- mx = info -> mx;
- my = info -> my;
-
- /* Define global and local grid parameters */
-
- dx = param->width/(PetscScalar) (mx); /* for periodic domains where the period is size width, then there are mx+1 grid points per period */
- dz = param->height/(PetscScalar) (my); /* ditto */
- dzdx2=dz*dz/dx/dx; /* grid aspect ratio squared */
- dz2=dz*dz;
-
-
- xints = info->xs; xinte = info->xs+info->xm;
- yints = info->ys; yinte = info->ys+info->ym;
-
- /* get time stepping parameters */
- t=tsCtx->t;
- dt=tsCtx->dt;
-
- /* if not full step, just calculate compaction rate for dt=0 */
- if ( tsCtx->takeFullStep) {
- halfDt=dt/2;
- } else {
- halfDt=0.;
- }
-
- /* zero the array pointer memory and unused k components*/
-
- PetscMemzero(col,5*sizeof(MatStencil));
- row.k = 0;
-
-
- /*
- porosity and compaction rate stencils for all points (periodic plus ghosts should require
- no boundary conditions)
- */
-
- for (j=yints; j<yinte; j++) {
- if (! (param->quiet)) {
- PetscPrintf(PETSC_COMM_WORLD,"row %d ",j);
- }
- for (i=xints; i<xinte; i++) {
-
- /* Calculate porosities and permeabilities at half nodes */
- phlfE = .5*(x[j][i].phi+x[j][i+1].phi);
- phlfW = .5*(x[j][i].phi+x[j][i-1].phi);
- phlfN = .5*(x[j][i].phi+x[j+1][i].phi);
- phlfS = .5*(x[j][i].phi+x[j-1][i].phi);
-
- khlfE = pow(phlfE,nExp);
- khlfW = pow(phlfW,nExp);
- khlfN = pow(phlfN,nExp);
- khlfS = pow(phlfS,nExp);
- sum = ( khlfN + khlfS + dzdx2*(khlfE + khlfW) + dz2);
-
- /*
- Jacobian entries for compaction rate equation (assuming permeabilities are passive)
- f[j][i].cmp = dzdx2*(khlfE*x[j][i+1].cmp + khlfW*x[j][i-1].cmp) +
- khlfN*x[j+1][i].cmp + khlfS*x[j-1][i].cmp + sum*x[j][i].cmp - dz*(khlfN-khlfS);
- */
- row.i = i; row.j = j; row.c = 0;
- v[0] = -khlfS; col[0].i = i; col[0].j = j-1; col[0].c = 0;
- v[1] = -dzdx2*khlfW; col[1].i = i-1; col[1].j = j; col[1].c = 0;
- v[2] = sum; col[2].i = i; col[2].j = j; col[2].c = 0;
- v[3] = -dzdx2*khlfE; col[3].i = i+1; col[3].j = j; col[3].c = 0;
- v[4] = -khlfN; col[4].i = i; col[4].j = j+1; col[4].c = 0;
- MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
-
- /*
- Jacobian entries for porosity update equation
- f[j][i].phi = x[j][i].phi - xold[j][i].phi - halfDt*(xold[j][i].cmp+x[j][i].cmp);
- */
- row.i = i; row.j = j; row.c = 1;
- v[0] = one; col[0].i = i; col[0].j = j; col[0].c = 1;
- v[1] = -halfDt; col[1].i = i; col[1].j = j; col[1].c = 0;
- MatSetValuesStencil(jac,1,&row,2,col,v,INSERT_VALUES);
}
}
@@ -575,29 +458,4 @@ PetscReal CubicInterp(PetscReal x, Petsc
return retval;
}
-/*---------------------------------------------------------------------*/
-#undef __FUNCT__
-#define __FUNCT__ "smoothRedBlack2D"
-PetscErrorCode smoothPorosityRedBlack2D(Field **x, int iMin, int iMax, int jMin, int jMax,int nIts)
-
-/*
- Performs nIts sweeps of red-black Gauss-Seidel to filter high frequency checkerboards
-
----------------------------------------------------------------------*/
-{
- PetscReal qtr=0.25;
- PetscErrorCode retval=0;
- int i, j, k, n;
-
-
- for (n=0; n< nIts; n++ ){
- for (k=0; k< 2; k++ ){
- for (j = jMin; j< jMax; j++ ) {
- for( i = iMin +k ; i < iMax; i+= 2 ) {
- x[j][i].phi = qtr*(x[j-1][i].phi + x[j][i-1].phi + x[j][i+1].phi + x[j+1][i].phi );
- }
- }
- }
- }
- return retval;
-}
+
diff -r b68f2b824135 -r 4adadc50962c src/Params.c
--- a/src/Params.c Sun Mar 16 22:44:05 2008 -0400
+++ b/src/Params.c Sun Mar 16 23:47:55 2008 -0400
@@ -52,12 +52,14 @@ int SetDomainAndWaveParams(PetscBag bag,
int ierr;
REG_INTG(bag,&p->mgLevels,5 ,"mglevels","total number of multi-grid levels, set with -mglevels");
REG_INTG(bag,&p->nDims,2 ,"nDims","<DO NOT SET> Grid dimensionality");
- REG_INTG(bag,&p->niCoarse,5 ,"ni_coarse","number of x points on coarsest level, set with -ni_coarse");
- REG_INTG(bag,&p->njCoarse,5 ,"nj_coarse","number of x points on coarsest level, set with -nj_coarse");
+ REG_INTG(bag,&p->niCoarse,4 ,"ni_coarse","number of x points on coarsest level, set with -ni_coarse");
+ REG_INTG(bag,&p->njCoarse,4 ,"nj_coarse","number of x points on coarsest level, set with -nj_coarse");
REG_INTG(bag,&p->xBoundaryFlag,0 ,"xBCFlag","Boundary condition in X, set with -xBCFlag (0: periodic, 1: reflection");
REG_INTG(bag,&p->zBoundaryFlag,0 ,"zBCFlag","Boundary condition in z, set with -zBCFlag (0: periodic, 1: phi=1,C=0 top, freeflux bottom");
REG_REAL(bag,&p->width,64.0 ,"width" ,"Width of domain, x-dir, set with -width");
REG_REAL(bag,&p->height,64.0 ,"height","height of domain, z-dir, set with -height");
+ REG_REAL(bag,&p->dx,1. ,"dx","<DO NOT SET> horizontal grid spacing (fine grid)");
+ REG_REAL(bag,&p->dz,1. ,"dz","<DO NOT SET> vertical grid spacing (fine grid)");
REG_REAL(bag,&p->nExp,3. ,"nExp","permeability exponent, set with -nExp ");
REG_REAL(bag,&p->amp,2. ,"amp","amplitude of initial solitary wave, with -amp ");
REG_REAL(bag,&p->z0,32.0 ,"z0","initial height of solitary wave, set with -z0 ");
diff -r b68f2b824135 -r 4adadc50962c src/PetscSimOutput.c
--- a/src/PetscSimOutput.c Sun Mar 16 22:44:05 2008 -0400
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,192 +0,0 @@
-/* ----------------------------------------------------------------------
- * PetscSimOutput.c
- *
- * Source file for PetscSimOutput library of functions to write
- * info files with matlab code for interpreting various PETSc
- * binary files.
- *
- * Note all "name" and "DAFieldName" variables must be Matlab Kosher-
- * i.e. no whitespace or illegal characters such as grouping
- * operators, quotations, math/boolean operators, etc.
- *
- * source: PetscSimOutput.c
- *
- * ----------------------------------------------------------------------*/
-
-
-
-#include "PetscSimOutput.h"
-
-
-/* ---------------------------------------------------------------------
- * PetscWriteOutputInitialize
- *
- * Open viewer and write matlab info file initialization.
- *
- * Input
- * ------------------------------
- * comm | mpi communicator
- * fname | filename
- * viewer | petsc viewer object
- *
- * ----------------------------------------------------------------------*/
-
-#undef __FUNCT__
-#define __FUNCT__ "PetscWriteOutputInitialize"
-PetscErrorCode PetscWriteOutputInitialize(MPI_Comm comm, const char *fname, PetscViewer *viewer)
-{
- PetscErrorCode ierr;
- FILE *info;
- PetscFunctionBegin;
- ierr = PetscViewerBinaryOpen(comm,fname,FILE_MODE_WRITE,viewer);CHKERRQ(ierr);
- ierr = PetscViewerBinaryGetInfoPointer(*viewer,&info);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscWriteOutputInitialize ---%\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ Set.filename = '%s';\n",fname);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ fd = fopen(Set.filename, 'r', 'ieee-be');\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ if (fd < 0) error('Cannot open %s, check for existence of file'); end\n",fname);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscWriteOutputInitialize ---%\n\n");CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-
-/* ---------------------------------------------------------------------
- * PetscWriteOutputFinalize
- *
- * Write matlab info file finalization and destroy viewer.
- *
- * Input
- * ------------------------------
- * viewer | petsc viewer object
- *
- * ----------------------------------------------------------------------*/
-
-#undef __FUNCT__
-#define __FUNCT__ "PetscWriteOutputFinalize"
-PetscErrorCode PetscWriteOutputFinalize(PetscViewer viewer)
-{
- PetscErrorCode ierr;
- FILE *info;
- MPI_Comm comm;
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
- ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscWriteOutputFinalize ---%\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ fclose(fd);\n");
- ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscWriteOutputFinalize ---%\n\n");CHKERRQ(ierr);
- ierr = PetscViewerFlush(viewer);CHKERRQ(ierr);
- ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-
-
-/* ---------------------------------------------------------------------
- * PetscWriteOutputBag
- *
- * Write matlab code to info file to read a PetscBag from binary.
- *
- * Input
- * ------------------------------
- * viewer | petsc viewer object
- * name | petsc bag name
- * bag | bag object of data to output
- *
- * ----------------------------------------------------------------------*/
-
-#undef __FUNCT__
-#define __FUNCT__ "PetscWriteOutputBag"
-PetscErrorCode PetscWriteOutputBag(PetscViewer viewer, const char *name, PetscBag bag)
-{
- PetscErrorCode ierr;
- FILE *info;
- MPI_Comm comm;
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
- ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
- ierr = PetscBagView(bag,viewer);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscWriteOutputBag ---%\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ Set.%s = PetscBinaryRead(fd);\n",name);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscWriteOutputBag ---%\n\n");CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-/* ---------------------------------------------------------------------
- * PetscWriteOutputVec
- *
- * Write matlab code to info file to read a (non DA) Vec from binary.
- *
- * Input
- * ------------------------------
- * viewer | petsc viewer object
- * name | name of variable to be written
- * vec | Vec of data to output
- *
- * ----------------------------------------------------------------------*/
-
-#undef __FUNCT__
-#define __FUNCT__ "PetscWriteOutputVec"
-PetscErrorCode PetscWriteOutputVec(PetscViewer viewer, const char* name, Vec vec)
-{
- PetscErrorCode ierr;
- FILE *info;
- MPI_Comm comm;
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
- ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
- ierr = VecView(vec,viewer);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscWriteOutputVec ---%\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ Set.%s = PetscBinaryRead(fd);\n",name);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscWriteOutputVec ---%\n\n");CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
-
-
-/* ---------------------------------------------------------------------
- * PetscWriteOutputVecDA
- *
- * Write matlab code to info file to read a DA's Vec from binary.
- * Note that this requires dof names have been set using DASetFieldName.
- *
- * Input
- * ------------------------------
- * viewer | petsc viewer object
- * name | names of variable to be written
- * vec | Vec of data to output
- * da | DA governing implementation of vec
- *
- * ----------------------------------------------------------------------*/
-
-#undef __FUNCT__
-#define __FUNCT__ "PetscWriteOutputVecDA"
-PetscErrorCode PetscWriteOutputVecDA(PetscViewer viewer, const char* name, Vec vec, DA da)
-{
- PetscErrorCode ierr;
- FILE *info;
- MPI_Comm comm;
- PetscInt dim,ni,nj,nk,pi,pj,pk,dof,n;
- PetscTruth flg;
- char *fieldname;
- PetscFunctionBegin;
- ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr);
- ierr = PetscViewerBinaryGetInfoPointer(viewer,&info);CHKERRQ(ierr);
- ierr = DAGetInfo(da,&dim,&ni,&nj,&nk,&pi,&pj,&pk,&dof,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
- ierr = VecView(vec,viewer);CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- begin code written by PetscWriteOutputVecDA ---%\n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%$$ tmp = PetscBinaryRead(fd); \n");CHKERRQ(ierr);
- if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d);\n",dof,ni);CHKERRQ(ierr); }
- if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d);\n",dof,ni,nj);CHKERRQ(ierr); }
- if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ tmp = reshape(tmp,%d,%d,%d,%d);\n",dof,ni,nj,nk);CHKERRQ(ierr); }
-
- for(n=0; n<dof; n++) {
- ierr = DAGetFieldName(da,n,&fieldname); CHKERRQ(ierr);
- ierr = PetscStrcmp(fieldname,"",&flg);CHKERRQ(ierr);
- if (!flg) {
- if (dim == 1) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
- if (dim == 2) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = squeeze(tmp(%d,:,:))';\n",name,fieldname,n+1);CHKERRQ(ierr); }
- if (dim == 3) { ierr = PetscFPrintf(comm,info,"%%$$ Set.%s.%s = permute(squeeze(tmp(%d,:,:,:)),[2 1 3]);\n",name,fieldname,n+1);CHKERRQ(ierr);}
- }
- }
- ierr = PetscFPrintf(comm,info,"%%$$ clear tmp; \n");CHKERRQ(ierr);
- ierr = PetscFPrintf(comm,info,"%%--- end code written by PetscWriteOutputVecDA ---%\n\n");CHKERRQ(ierr);
- PetscFunctionReturn(0);
-}
diff -r b68f2b824135 -r 4adadc50962c src/SolitaryWave2DSL.c
--- a/src/SolitaryWave2DSL.c Sun Mar 16 22:44:05 2008 -0400
+++ b/src/SolitaryWave2DSL.c Sun Mar 16 23:47:55 2008 -0400
@@ -108,16 +108,18 @@ int main(int argc,char **argv)
ierr = DAGetInfo(DMMGGetDA(dmmg),0,&grid.niFine,&grid.njFine,0,0,0,0,0,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
/* set correct grid spacing depending on periodicity of boundary conditions */
if ( param->xBoundaryFlag == 0 ) {
- grid.dx = param->width / (PetscReal) grid.niFine;
+ param->dx = param->width / (PetscReal) grid.niFine;
} else {
- grid.dx = param->width / (PetscReal) (grid.niFine - 1);
+ param->dx = param->width / (PetscReal) (grid.niFine - 1);
}
if ( param->zBoundaryFlag == 0 ) {
- grid.dz = param->height / (PetscReal) grid.njFine;
+ param->dz = param->height / (PetscReal) grid.njFine;
} else {
- grid.dz = param->height / (PetscReal) (grid.njFine -1);
- }
- tsCtx.W0Scaled = param->W0/grid.dz;
+ param->dz = param->height / (PetscReal) (grid.njFine -1);
+ }
+ grid.dx = param->dx;
+ grid.dz = param->dz;
+ tsCtx.W0Scaled = param->W0/param->dz;
/* create DA for characteristics solve (stencil width =2) */
ierr = DACreate2d(comm,grid.periodic,grid.stencilC,grid.niFine,grid.njFine,PETSC_DECIDE,PETSC_DECIDE,
@@ -150,11 +152,7 @@ int main(int argc,char **argv)
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Set up the SNES solver with callback functions (including local Jacobian)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
- if (param->useApproxJacobian) {
- ierr = DMMGSetSNESLocal(dmmg,FormFunctionLocal,FormJacobianLocalApprox,0,0);CHKERRQ(ierr);
- } else {
- ierr = DMMGSetSNESLocal(dmmg,FormFunctionLocal,FormJacobianLocal,0,0);CHKERRQ(ierr);
- }
+ ierr = DMMGSetSNESLocal(dmmg,FormFunctionLocal,FormJacobianLocal,0,0);CHKERRQ(ierr);
ierr = DMMGSetInitialGuess(dmmg,FormInitialGuess);CHKERRQ(ierr);
snes = DMMGGetSNES(dmmg);CHKERRQ(ierr);
diff -r b68f2b824135 -r 4adadc50962c src/SolitaryWave2DSL.h
--- a/src/SolitaryWave2DSL.h Sun Mar 16 22:44:05 2008 -0400
+++ b/src/SolitaryWave2DSL.h Sun Mar 16 23:47:55 2008 -0400
@@ -49,6 +49,7 @@ typedef struct {
int mgLevels, nDims; /* multi-grid parameters */
int niCoarse, njCoarse; /* coarse grid parameters */
PetscReal width, height; /* width and height of box in compaction lengths */
+ PetscReal dx, dz; /* horizontal and vertical grid-spacing */
PetscTruth useFieldSplitPC; /* Boolean to initialize fieldsplitPC...or not */
PetscTruth useApproxJacobian; /* Boolean to toggle full or approximate Jacobian */
/* Boundary Condtion Flags */
@@ -102,7 +103,6 @@ PetscErrorCode FormInitialGuess(DMMG,Vec
PetscErrorCode FormInitialGuess(DMMG,Vec);
PetscErrorCode FormFunctionLocal(DALocalInfo*,Field**,Field**,void*);
PetscErrorCode FormJacobianLocal(DALocalInfo*,Field**,Mat,void*);
-PetscErrorCode FormJacobianLocalApprox(DALocalInfo*,Field**,Mat,void*);
/* Utility routines */
int SetParams(AppCtx*);
diff -r b68f2b824135 -r 4adadc50962c src/makeOutput.c
--- a/src/makeOutput.c Sun Mar 16 22:44:05 2008 -0400
+++ b/src/makeOutput.c Sun Mar 16 23:47:55 2008 -0400
@@ -9,7 +9,8 @@ int DoOutput(DMMG *dmmg, int nInt, int n
AppCtx *user = (AppCtx*)dmmg[0]->user;
PetscBag bag = user->bag;
Parameter *param;
- TStepCtx *tsCtx = user->tsCtx;
+ TStepCtx *tsCtx = user->tsCtx;
+
DA da = DMMGGetDA(dmmg);
int ierr;
More information about the cig-commits
mailing list