Fixed a problem with Set(Box)Values when called after AddTo(Box)Values.

Set should zero out entries in ghost zones to ensure that the accumulation
that happens at assembly time is correct.
This commit is contained in:
falgout 2008-05-05 19:32:49 +00:00
parent 743b325637
commit e512671109
5 changed files with 430 additions and 17 deletions

View File

@ -290,6 +290,8 @@ hypre_SStructPMatrixSetValues( hypre_SStructPMatrix *pmatrix,
int *smap = hypre_SStructPMatrixSMap(pmatrix, var);
int *vars = hypre_SStructStencilVars(stencil);
hypre_StructMatrix *smatrix;
hypre_BoxArray *grid_boxes;
hypre_Box *box;
int *sentries;
int i;
@ -305,13 +307,12 @@ hypre_SStructPMatrixSetValues( hypre_SStructPMatrix *pmatrix,
hypre_StructMatrixSetValues(smatrix, index, nentries, sentries, values,
action, -1, 0);
/* set values outside the grid in ghost zones (for AddTo and Get) */
/* set (AddTo/Get) or clear (Set) values outside the grid in ghost zones */
if (action != 0)
{
/* AddTo/Get */
hypre_SStructPGrid *pgrid = hypre_SStructPMatrixPGrid(pmatrix);
hypre_Index varoffset;
hypre_BoxArray *grid_boxes;
hypre_Box *box;
int done = 0;
grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(smatrix));
@ -358,6 +359,25 @@ hypre_SStructPMatrixSetValues( hypre_SStructPMatrix *pmatrix,
}
}
}
else
{
/* Set */
grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(smatrix));
hypre_ForBoxI(i, grid_boxes)
{
box = hypre_BoxArrayBox(grid_boxes, i);
if ((hypre_IndexX(index) < hypre_BoxIMinX(box)) ||
(hypre_IndexX(index) > hypre_BoxIMaxX(box)) ||
(hypre_IndexY(index) < hypre_BoxIMinY(box)) ||
(hypre_IndexY(index) > hypre_BoxIMaxY(box)) ||
(hypre_IndexZ(index) < hypre_BoxIMinZ(box)) ||
(hypre_IndexZ(index) > hypre_BoxIMaxZ(box)) )
{
hypre_StructMatrixClearValues(smatrix, index, nentries, sentries, i, 1);
}
}
}
return hypre_error_flag;
}
@ -383,10 +403,11 @@ hypre_SStructPMatrixSetBoxValues( hypre_SStructPMatrix *pmatrix,
int *smap = hypre_SStructPMatrixSMap(pmatrix, var);
int *vars = hypre_SStructStencilVars(stencil);
hypre_StructMatrix *smatrix;
hypre_BoxArray *grid_boxes;
hypre_Box *box;
hypre_Box *value_box;
int *sentries;
int i;
int i, j;
smatrix = hypre_SStructPMatrixSMatrix(pmatrix, var, vars[entries[0]]);
@ -405,15 +426,14 @@ hypre_SStructPMatrixSetBoxValues( hypre_SStructPMatrix *pmatrix,
hypre_StructMatrixSetBoxValues(smatrix, box, value_box, nentries, sentries,
values, action, -1, 0);
/* set values outside the grid in ghost zones (for AddTo and Get) */
/* set (AddTo/Get) or clear (Set) values outside the grid in ghost zones */
if (action != 0)
{
/* AddTo/Get */
hypre_SStructPGrid *pgrid = hypre_SStructPMatrixPGrid(pmatrix);
hypre_Index varoffset;
hypre_BoxArray *grid_boxes;
hypre_BoxArray *left_boxes, *done_boxes, *temp_boxes;
hypre_Box *left_box, *done_box, *int_box;
int j;
hypre_SStructVariableGetOffset(hypre_SStructPGridVarType(pgrid, var),
hypre_SStructPGridNDim(pgrid), varoffset);
@ -458,6 +478,30 @@ hypre_SStructPMatrixSetBoxValues( hypre_SStructPMatrix *pmatrix,
hypre_BoxArrayDestroy(done_boxes);
hypre_BoxArrayDestroy(temp_boxes);
}
else
{
/* Set */
hypre_BoxArray *diff_boxes;
hypre_Box *grid_box, *diff_box;
grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(smatrix));
diff_boxes = hypre_BoxArrayCreate(0);
hypre_ForBoxI(i, grid_boxes)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
hypre_BoxArraySetSize(diff_boxes, 0);
hypre_SubtractBoxes(box, grid_box, diff_boxes);
hypre_ForBoxI(j, diff_boxes)
{
diff_box = hypre_BoxArrayBox(diff_boxes, j);
hypre_StructMatrixClearBoxValues(smatrix, diff_box, nentries, sentries,
i, 1);
}
}
hypre_BoxArrayDestroy(diff_boxes);
}
hypre_BoxDestroy(box);

View File

@ -179,19 +179,20 @@ hypre_SStructPVectorSetValues( hypre_SStructPVector *pvector,
int action )
{
hypre_StructVector *svector = hypre_SStructPVectorSVector(pvector, var);
hypre_BoxArray *grid_boxes;
hypre_Box *box;
int i;
/* set values inside the grid */
hypre_StructVectorSetValues(svector, index, value, action, -1, 0);
/* set values outside the grid in ghost zones (for AddTo and Get) */
/* set (AddTo/Get) or clear (Set) values outside the grid in ghost zones */
if (action != 0)
{
hypre_SStructPGrid *pgrid = hypre_SStructPVectorPGrid(pvector);
hypre_Index varoffset;
hypre_BoxArray *grid_boxes;
hypre_Box *box;
int i;
int done = 0;
/* AddTo/Get */
hypre_SStructPGrid *pgrid = hypre_SStructPVectorPGrid(pvector);
hypre_Index varoffset;
int done = 0;
grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(svector));
@ -236,6 +237,25 @@ hypre_SStructPVectorSetValues( hypre_SStructPVector *pvector,
}
}
}
else
{
/* Set */
grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(svector));
hypre_ForBoxI(i, grid_boxes)
{
box = hypre_BoxArrayBox(grid_boxes, i);
if ((hypre_IndexX(index) < hypre_BoxIMinX(box)) ||
(hypre_IndexX(index) > hypre_BoxIMaxX(box)) ||
(hypre_IndexY(index) < hypre_BoxIMinY(box)) ||
(hypre_IndexY(index) > hypre_BoxIMaxY(box)) ||
(hypre_IndexZ(index) < hypre_BoxIMinZ(box)) ||
(hypre_IndexZ(index) > hypre_BoxIMaxZ(box)) )
{
hypre_StructVectorClearValues(svector, index, i, 1);
}
}
}
return hypre_error_flag;
}
@ -255,8 +275,10 @@ hypre_SStructPVectorSetBoxValues( hypre_SStructPVector *pvector,
int action )
{
hypre_StructVector *svector = hypre_SStructPVectorSVector(pvector, var);
hypre_BoxArray *grid_boxes;
hypre_Box *box;
hypre_Box *value_box;
int i, j;
box = hypre_BoxCreate();
hypre_CopyIndex(ilower, hypre_BoxIMin(box));
@ -266,15 +288,14 @@ hypre_SStructPVectorSetBoxValues( hypre_SStructPVector *pvector,
/* set values inside the grid */
hypre_StructVectorSetBoxValues(svector, box, value_box, values, action, -1, 0);
/* set values outside the grid in ghost zones (for AddTo and Get) */
/* set (AddTo/Get) or clear (Set) values outside the grid in ghost zones */
if (action != 0)
{
/* AddTo/Get */
hypre_SStructPGrid *pgrid = hypre_SStructPVectorPGrid(pvector);
hypre_Index varoffset;
hypre_BoxArray *grid_boxes;
hypre_BoxArray *left_boxes, *done_boxes, *temp_boxes;
hypre_Box *left_box, *done_box, *int_box;
int i, j;
hypre_SStructVariableGetOffset(hypre_SStructPGridVarType(pgrid, var),
hypre_SStructPGridNDim(pgrid), varoffset);
@ -318,6 +339,29 @@ hypre_SStructPVectorSetBoxValues( hypre_SStructPVector *pvector,
hypre_BoxArrayDestroy(done_boxes);
hypre_BoxArrayDestroy(temp_boxes);
}
else
{
/* Set */
hypre_BoxArray *diff_boxes;
hypre_Box *grid_box, *diff_box;
grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(svector));
diff_boxes = hypre_BoxArrayCreate(0);
hypre_ForBoxI(i, grid_boxes)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
hypre_BoxArraySetSize(diff_boxes, 0);
hypre_SubtractBoxes(box, grid_box, diff_boxes);
hypre_ForBoxI(j, diff_boxes)
{
diff_box = hypre_BoxArrayBox(diff_boxes, j);
hypre_StructVectorClearBoxValues(svector, diff_box, i, 1);
}
}
hypre_BoxArrayDestroy(diff_boxes);
}
hypre_BoxDestroy(box);

View File

@ -2349,6 +2349,8 @@ int hypre_StructMatrixInitialize ( hypre_StructMatrix *matrix );
int hypre_StructMatrixSetValues ( hypre_StructMatrix *matrix , hypre_Index grid_index , int num_stencil_indices , int *stencil_indices , double *values , int action , int boxnum , int outside );
int hypre_StructMatrixSetBoxValues ( hypre_StructMatrix *matrix , hypre_Box *set_box , hypre_Box *value_box , int num_stencil_indices , int *stencil_indices , double *values , int action , int boxnum , int outside );
int hypre_StructMatrixSetConstantValues ( hypre_StructMatrix *matrix , int num_stencil_indices , int *stencil_indices , double *values , int action );
int hypre_StructMatrixClearValues ( hypre_StructMatrix *matrix , hypre_Index grid_index , int num_stencil_indices , int *stencil_indices , int boxnum , int outside );
int hypre_StructMatrixClearBoxValues ( hypre_StructMatrix *matrix , hypre_Box *clear_box , int num_stencil_indices , int *stencil_indices , int boxnum , int outside );
int hypre_StructMatrixAssemble ( hypre_StructMatrix *matrix );
int hypre_StructMatrixSetNumGhost ( hypre_StructMatrix *matrix , int *num_ghost );
int hypre_StructMatrixSetConstantCoefficient ( hypre_StructMatrix *matrix , int constant_coefficient );
@ -2393,6 +2395,8 @@ int hypre_StructVectorInitializeData ( hypre_StructVector *vector , double *data
int hypre_StructVectorInitialize ( hypre_StructVector *vector );
int hypre_StructVectorSetValues ( hypre_StructVector *vector , hypre_Index grid_index , double *values , int action , int boxnum , int outside );
int hypre_StructVectorSetBoxValues ( hypre_StructVector *vector , hypre_Box *set_box , hypre_Box *value_box , double *values , int action , int boxnum , int outside );
int hypre_StructVectorClearValues ( hypre_StructVector *vector , hypre_Index grid_index , int boxnum , int outside );
int hypre_StructVectorClearBoxValues ( hypre_StructVector *vector , hypre_Box *clear_box , int boxnum , int outside );
int hypre_StructVectorSetNumGhost ( hypre_StructVector *vector , int *num_ghost );
int hypre_StructVectorAssemble ( hypre_StructVector *vector );
int hypre_StructVectorCopy ( hypre_StructVector *x , hypre_StructVector *y );

View File

@ -1057,6 +1057,177 @@ hypre_StructMatrixSetConstantValues( hypre_StructMatrix *matrix,
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* (outside > 0): clear values possibly outside of the grid extents
* (outside = 0): clear values only inside the grid extents
*--------------------------------------------------------------------------*/
int
hypre_StructMatrixClearValues( hypre_StructMatrix *matrix,
hypre_Index grid_index,
int num_stencil_indices,
int *stencil_indices,
int boxnum,
int outside )
{
hypre_BoxArray *grid_boxes;
hypre_Box *grid_box;
double *matp;
int i, s, istart, istop;
/*-----------------------------------------------------------------------
* Initialize some things
*-----------------------------------------------------------------------*/
if (outside > 0)
{
grid_boxes = hypre_StructMatrixDataSpace(matrix);
}
else
{
grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
}
if (boxnum < 0)
{
istart = 0;
istop = hypre_BoxArraySize(grid_boxes);
}
else
{
istart = boxnum;
istop = istart + 1;
}
/*-----------------------------------------------------------------------
* Clear the matrix coefficients
*-----------------------------------------------------------------------*/
for (i = istart; i < istop; i++)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(grid_box)) &&
(hypre_IndexX(grid_index) <= hypre_BoxIMaxX(grid_box)) &&
(hypre_IndexY(grid_index) >= hypre_BoxIMinY(grid_box)) &&
(hypre_IndexY(grid_index) <= hypre_BoxIMaxY(grid_box)) &&
(hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(grid_box)) &&
(hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(grid_box)) )
{
for (s = 0; s < num_stencil_indices; s++)
{
matp = hypre_StructMatrixBoxDataValue(matrix, i,
stencil_indices[s],
grid_index);
*matp = 0.0;
}
}
}
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* (outside > 0): clear values possibly outside of the grid extents
* (outside = 0): clear values only inside the grid extents
*--------------------------------------------------------------------------*/
int
hypre_StructMatrixClearBoxValues( hypre_StructMatrix *matrix,
hypre_Box *clear_box,
int num_stencil_indices,
int *stencil_indices,
int boxnum,
int outside )
{
hypre_BoxArray *grid_boxes;
hypre_Box *grid_box;
hypre_Box *int_box;
hypre_BoxArray *data_space;
hypre_Box *data_box;
hypre_IndexRef data_start;
hypre_Index data_stride;
int datai;
double *datap;
hypre_Index loop_size;
int i, s, istart, istop;
int loopi, loopj, loopk;
/*-----------------------------------------------------------------------
* Initialize some things
*-----------------------------------------------------------------------*/
if (outside > 0)
{
grid_boxes = hypre_StructMatrixDataSpace(matrix);
}
else
{
grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
}
data_space = hypre_StructMatrixDataSpace(matrix);
if (boxnum < 0)
{
istart = 0;
istop = hypre_BoxArraySize(grid_boxes);
}
else
{
istart = boxnum;
istop = istart + 1;
}
/*-----------------------------------------------------------------------
* Clear the matrix coefficients
*-----------------------------------------------------------------------*/
hypre_SetIndex(data_stride, 1, 1, 1);
int_box = hypre_BoxCreate();
for (i = istart; i < istop; i++)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
data_box = hypre_BoxArrayBox(data_space, i);
hypre_IntersectBoxes(clear_box, grid_box, int_box);
/* if there was an intersection */
if (int_box)
{
data_start = hypre_BoxIMin(int_box);
for (s = 0; s < num_stencil_indices; s++)
{
datap = hypre_StructMatrixBoxData(matrix, i,
stencil_indices[s]);
hypre_BoxGetSize(int_box, loop_size);
hypre_BoxLoop1Begin(loop_size,
data_box,data_start,data_stride,datai);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,datai
#include "hypre_box_smp_forloop.h"
hypre_BoxLoop1For(loopi, loopj, loopk, datai)
{
datap[datai] = 0.0;
}
hypre_BoxLoop1End(datai);
}
}
}
hypre_BoxDestroy(int_box);
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
*--------------------------------------------------------------------------*/

View File

@ -426,6 +426,156 @@ hypre_StructVectorSetBoxValues( hypre_StructVector *vector,
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* (outside > 0): clear values possibly outside of the grid extents
* (outside = 0): clear values only inside the grid extents
*--------------------------------------------------------------------------*/
int
hypre_StructVectorClearValues( hypre_StructVector *vector,
hypre_Index grid_index,
int boxnum,
int outside )
{
hypre_BoxArray *grid_boxes;
hypre_Box *grid_box;
double *vecp;
int i, istart, istop;
if (outside > 0)
{
grid_boxes = hypre_StructVectorDataSpace(vector);
}
else
{
grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
}
if (boxnum < 0)
{
istart = 0;
istop = hypre_BoxArraySize(grid_boxes);
}
else
{
istart = boxnum;
istop = istart + 1;
}
for (i = istart; i < istop; i++)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(grid_box)) &&
(hypre_IndexX(grid_index) <= hypre_BoxIMaxX(grid_box)) &&
(hypre_IndexY(grid_index) >= hypre_BoxIMinY(grid_box)) &&
(hypre_IndexY(grid_index) <= hypre_BoxIMaxY(grid_box)) &&
(hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(grid_box)) &&
(hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(grid_box)) )
{
vecp = hypre_StructVectorBoxDataValue(vector, i, grid_index);
*vecp = 0.0;
}
}
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* (outside > 0): clear values possibly outside of the grid extents
* (outside = 0): clear values only inside the grid extents
*--------------------------------------------------------------------------*/
int
hypre_StructVectorClearBoxValues( hypre_StructVector *vector,
hypre_Box *clear_box,
int boxnum,
int outside )
{
hypre_BoxArray *grid_boxes;
hypre_Box *grid_box;
hypre_Box *int_box;
hypre_BoxArray *data_space;
hypre_Box *data_box;
hypre_IndexRef data_start;
hypre_Index data_stride;
int datai;
double *datap;
hypre_Index loop_size;
int i, istart, istop;
int loopi, loopj, loopk;
/*-----------------------------------------------------------------------
* Initialize some things
*-----------------------------------------------------------------------*/
if (outside > 0)
{
grid_boxes = hypre_StructVectorDataSpace(vector);
}
else
{
grid_boxes = hypre_StructGridBoxes(hypre_StructVectorGrid(vector));
}
data_space = hypre_StructVectorDataSpace(vector);
if (boxnum < 0)
{
istart = 0;
istop = hypre_BoxArraySize(grid_boxes);
}
else
{
istart = boxnum;
istop = istart + 1;
}
/*-----------------------------------------------------------------------
* Set the vector coefficients
*-----------------------------------------------------------------------*/
hypre_SetIndex(data_stride, 1, 1, 1);
int_box = hypre_BoxCreate();
for (i = istart; i < istop; i++)
{
grid_box = hypre_BoxArrayBox(grid_boxes, i);
data_box = hypre_BoxArrayBox(data_space, i);
hypre_IntersectBoxes(clear_box, grid_box, int_box);
/* if there was an intersection */
if (int_box)
{
data_start = hypre_BoxIMin(int_box);
datap = hypre_StructVectorBoxData(vector, i);
hypre_BoxGetSize(int_box, loop_size);
hypre_BoxLoop1Begin(loop_size,
data_box,data_start,data_stride,datai);
#define HYPRE_BOX_SMP_PRIVATE loopk,loopi,loopj,datai
#include "hypre_box_smp_forloop.h"
hypre_BoxLoop1For(loopi, loopj, loopk, datai)
{
datap[datai] = 0.0;
}
hypre_BoxLoop1End(datai);
}
}
hypre_BoxDestroy(int_box);
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
*--------------------------------------------------------------------------*/