Adding AMS support for 1D domains (#269)
This commit is contained in:
Mark L. Stowell
GitHub
parent
5f8472b05c
commit
970f75821f
@ -763,7 +763,7 @@ HYPRE_Int hypre_AMSSetDimension(void *solver,
|
||||
{
|
||||
hypre_AMSData *ams_data = (hypre_AMSData *) solver;
|
||||
|
||||
if (dim != 2 && dim != 3)
|
||||
if (dim != 1 && dim != 2 && dim != 3)
|
||||
hypre_error_in_arg(2);
|
||||
|
||||
ams_data -> dim = dim;
|
||||
@ -1213,7 +1213,8 @@ HYPRE_Int hypre_AMSComputePi(hypre_ParCSRMatrix *A,
|
||||
hypre_ParCSRMatrixInitialize(Pi);
|
||||
|
||||
Gx_data = hypre_VectorData(hypre_ParVectorLocalVector(Gx));
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim >= 2)
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim == 3)
|
||||
Gz_data = hypre_VectorData(hypre_ParVectorLocalVector(Gz));
|
||||
|
||||
@ -1243,7 +1244,8 @@ HYPRE_Int hypre_AMSComputePi(hypre_ParCSRMatrix *A,
|
||||
for (j = G_diag_I[i]; j < G_diag_I[i+1]; j++)
|
||||
{
|
||||
*Pi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gx_data[i];
|
||||
*Pi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim >= 2)
|
||||
*Pi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim == 3)
|
||||
*Pi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gz_data[i];
|
||||
}
|
||||
@ -1280,7 +1282,8 @@ HYPRE_Int hypre_AMSComputePi(hypre_ParCSRMatrix *A,
|
||||
for (j = G_offd_I[i]; j < G_offd_I[i+1]; j++)
|
||||
{
|
||||
*Pi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gx_data[i];
|
||||
*Pi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim >= 2)
|
||||
*Pi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim == 3)
|
||||
*Pi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gz_data[i];
|
||||
}
|
||||
@ -1349,18 +1352,21 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
hypre_ParCSRMatrixOwnsColStarts(Pix) = 0;
|
||||
hypre_ParCSRMatrixInitialize(Pix);
|
||||
|
||||
Piy = hypre_ParCSRMatrixCreate(comm,
|
||||
global_num_rows,
|
||||
global_num_cols,
|
||||
row_starts,
|
||||
col_starts,
|
||||
num_cols_offd,
|
||||
num_nonzeros_diag,
|
||||
num_nonzeros_offd);
|
||||
hypre_ParCSRMatrixOwnsData(Piy) = 1;
|
||||
hypre_ParCSRMatrixOwnsRowStarts(Piy) = 0;
|
||||
hypre_ParCSRMatrixOwnsColStarts(Piy) = 0;
|
||||
hypre_ParCSRMatrixInitialize(Piy);
|
||||
if (dim >= 2)
|
||||
{
|
||||
Piy = hypre_ParCSRMatrixCreate(comm,
|
||||
global_num_rows,
|
||||
global_num_cols,
|
||||
row_starts,
|
||||
col_starts,
|
||||
num_cols_offd,
|
||||
num_nonzeros_diag,
|
||||
num_nonzeros_offd);
|
||||
hypre_ParCSRMatrixOwnsData(Piy) = 1;
|
||||
hypre_ParCSRMatrixOwnsRowStarts(Piy) = 0;
|
||||
hypre_ParCSRMatrixOwnsColStarts(Piy) = 0;
|
||||
hypre_ParCSRMatrixInitialize(Piy);
|
||||
}
|
||||
|
||||
if (dim == 3)
|
||||
{
|
||||
@ -1379,7 +1385,8 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
}
|
||||
|
||||
Gx_data = hypre_VectorData(hypre_ParVectorLocalVector(Gx));
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim >= 2)
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim == 3)
|
||||
Gz_data = hypre_VectorData(hypre_ParVectorLocalVector(Gz));
|
||||
|
||||
@ -1431,7 +1438,7 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
*Piz_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gz_data[i];
|
||||
}
|
||||
}
|
||||
else
|
||||
else if (dim == 2)
|
||||
{
|
||||
hypre_CSRMatrix *G_diag = hypre_ParCSRMatrixDiag(G);
|
||||
HYPRE_Int *G_diag_I = hypre_CSRMatrixI(G_diag);
|
||||
@ -1470,6 +1477,37 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
*Piy_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gy_data[i];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
hypre_CSRMatrix *G_diag = hypre_ParCSRMatrixDiag(G);
|
||||
HYPRE_Int *G_diag_I = hypre_CSRMatrixI(G_diag);
|
||||
HYPRE_Int *G_diag_J = hypre_CSRMatrixJ(G_diag);
|
||||
HYPRE_Real *G_diag_data = hypre_CSRMatrixData(G_diag);
|
||||
|
||||
HYPRE_Int G_diag_nrows = hypre_CSRMatrixNumRows(G_diag);
|
||||
HYPRE_Int G_diag_nnz = hypre_CSRMatrixNumNonzeros(G_diag);
|
||||
|
||||
hypre_CSRMatrix *Pix_diag = hypre_ParCSRMatrixDiag(Pix);
|
||||
HYPRE_Int *Pix_diag_I = hypre_CSRMatrixI(Pix_diag);
|
||||
HYPRE_Int *Pix_diag_J = hypre_CSRMatrixJ(Pix_diag);
|
||||
HYPRE_Real *Pix_diag_data = hypre_CSRMatrixData(Pix_diag);
|
||||
|
||||
for (i = 0; i < G_diag_nrows+1; i++)
|
||||
{
|
||||
Pix_diag_I[i] = G_diag_I[i];
|
||||
}
|
||||
|
||||
for (i = 0; i < G_diag_nnz; i++)
|
||||
{
|
||||
Pix_diag_J[i] = G_diag_J[i];
|
||||
}
|
||||
|
||||
for (i = 0; i < G_diag_nrows; i++)
|
||||
for (j = G_diag_I[i]; j < G_diag_I[i+1]; j++)
|
||||
{
|
||||
*Pix_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gx_data[i];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* Fill-in the off-diagonal part */
|
||||
@ -1534,7 +1572,7 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
Piz_cmap[i] = G_cmap[i];
|
||||
}
|
||||
}
|
||||
else
|
||||
else if (dim == 2)
|
||||
{
|
||||
hypre_CSRMatrix *G_offd = hypre_ParCSRMatrixOffd(G);
|
||||
HYPRE_Int *G_offd_I = hypre_CSRMatrixI(G_offd);
|
||||
@ -1585,10 +1623,52 @@ HYPRE_Int hypre_AMSComputePixyz(hypre_ParCSRMatrix *A,
|
||||
Piy_cmap[i] = G_cmap[i];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
hypre_CSRMatrix *G_offd = hypre_ParCSRMatrixOffd(G);
|
||||
HYPRE_Int *G_offd_I = hypre_CSRMatrixI(G_offd);
|
||||
HYPRE_Int *G_offd_J = hypre_CSRMatrixJ(G_offd);
|
||||
HYPRE_Real *G_offd_data = hypre_CSRMatrixData(G_offd);
|
||||
|
||||
HYPRE_Int G_offd_nrows = hypre_CSRMatrixNumRows(G_offd);
|
||||
HYPRE_Int G_offd_ncols = hypre_CSRMatrixNumCols(G_offd);
|
||||
HYPRE_Int G_offd_nnz = hypre_CSRMatrixNumNonzeros(G_offd);
|
||||
|
||||
hypre_CSRMatrix *Pix_offd = hypre_ParCSRMatrixOffd(Pix);
|
||||
HYPRE_Int *Pix_offd_I = hypre_CSRMatrixI(Pix_offd);
|
||||
HYPRE_Int *Pix_offd_J = hypre_CSRMatrixJ(Pix_offd);
|
||||
HYPRE_Real *Pix_offd_data = hypre_CSRMatrixData(Pix_offd);
|
||||
|
||||
HYPRE_BigInt *G_cmap = hypre_ParCSRMatrixColMapOffd(G);
|
||||
HYPRE_BigInt *Pix_cmap = hypre_ParCSRMatrixColMapOffd(Pix);
|
||||
|
||||
if (G_offd_ncols)
|
||||
for (i = 0; i < G_offd_nrows+1; i++)
|
||||
{
|
||||
Pix_offd_I[i] = G_offd_I[i];
|
||||
}
|
||||
|
||||
for (i = 0; i < G_offd_nnz; i++)
|
||||
{
|
||||
Pix_offd_J[i] = G_offd_J[i];
|
||||
}
|
||||
|
||||
for (i = 0; i < G_offd_nrows; i++)
|
||||
for (j = G_offd_I[i]; j < G_offd_I[i+1]; j++)
|
||||
{
|
||||
*Pix_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gx_data[i];
|
||||
}
|
||||
|
||||
for (i = 0; i < G_offd_ncols; i++)
|
||||
{
|
||||
Pix_cmap[i] = G_cmap[i];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
*Pix_ptr = Pix;
|
||||
*Piy_ptr = Piy;
|
||||
if (dim >= 2)
|
||||
*Piy_ptr = Piy;
|
||||
if (dim == 3)
|
||||
*Piz_ptr = Piz;
|
||||
|
||||
@ -1653,7 +1733,8 @@ HYPRE_Int hypre_AMSComputeGPi(hypre_ParCSRMatrix *A,
|
||||
hypre_ParCSRMatrixInitialize(GPi);
|
||||
|
||||
Gx_data = hypre_VectorData(hypre_ParVectorLocalVector(Gx));
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim >= 3)
|
||||
Gy_data = hypre_VectorData(hypre_ParVectorLocalVector(Gy));
|
||||
if (dim == 4)
|
||||
Gz_data = hypre_VectorData(hypre_ParVectorLocalVector(Gz));
|
||||
|
||||
@ -1684,7 +1765,8 @@ HYPRE_Int hypre_AMSComputeGPi(hypre_ParCSRMatrix *A,
|
||||
{
|
||||
*GPi_diag_data++ = G_diag_data[j];
|
||||
*GPi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gx_data[i];
|
||||
*GPi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim >= 3)
|
||||
*GPi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim == 4)
|
||||
*GPi_diag_data++ = fabs(G_diag_data[j]) * 0.5 * Gz_data[i];
|
||||
}
|
||||
@ -1722,7 +1804,8 @@ HYPRE_Int hypre_AMSComputeGPi(hypre_ParCSRMatrix *A,
|
||||
{
|
||||
*GPi_offd_data++ = G_offd_data[j];
|
||||
*GPi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gx_data[i];
|
||||
*GPi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim >= 3)
|
||||
*GPi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gy_data[i];
|
||||
if (dim == 4)
|
||||
*GPi_offd_data++ = fabs(G_offd_data[j]) * 0.5 * Gz_data[i];
|
||||
}
|
||||
@ -1946,20 +2029,25 @@ HYPRE_Int hypre_AMSSetup(void *solver,
|
||||
|
||||
/* If not given, compute Gx, Gy and Gz */
|
||||
{
|
||||
if (ams_data -> x != NULL && ams_data -> y != NULL &&
|
||||
(ams_data -> dim == 2 || ams_data -> z != NULL))
|
||||
if (ams_data -> x != NULL &&
|
||||
(ams_data -> dim == 1 || ams_data -> y != NULL) &&
|
||||
(ams_data -> dim <= 2 || ams_data -> z != NULL))
|
||||
input_info = 1;
|
||||
|
||||
if (ams_data -> Gx != NULL && ams_data -> Gy != NULL &&
|
||||
(ams_data -> dim == 2 || ams_data -> Gz != NULL))
|
||||
if (ams_data -> Gx != NULL &&
|
||||
(ams_data -> dim == 1 || ams_data -> Gy != NULL) &&
|
||||
(ams_data -> dim <= 2 || ams_data -> Gz != NULL))
|
||||
input_info = 2;
|
||||
|
||||
if (input_info == 1)
|
||||
{
|
||||
ams_data -> Gx = hypre_ParVectorInRangeOf(ams_data -> G);
|
||||
hypre_ParCSRMatrixMatvec (1.0, ams_data -> G, ams_data -> x, 0.0, ams_data -> Gx);
|
||||
ams_data -> Gy = hypre_ParVectorInRangeOf(ams_data -> G);
|
||||
hypre_ParCSRMatrixMatvec (1.0, ams_data -> G, ams_data -> y, 0.0, ams_data -> Gy);
|
||||
if (ams_data -> dim >= 2)
|
||||
{
|
||||
ams_data -> Gy = hypre_ParVectorInRangeOf(ams_data -> G);
|
||||
hypre_ParCSRMatrixMatvec (1.0, ams_data -> G, ams_data -> y, 0.0, ams_data -> Gy);
|
||||
}
|
||||
if (ams_data -> dim == 3)
|
||||
{
|
||||
ams_data -> Gz = hypre_ParVectorInRangeOf(ams_data -> G);
|
||||
@ -2006,7 +2094,8 @@ HYPRE_Int hypre_AMSSetup(void *solver,
|
||||
if (input_info == 1 && ams_data -> cycle_type != 9)
|
||||
{
|
||||
hypre_ParVectorDestroy(ams_data -> Gx);
|
||||
hypre_ParVectorDestroy(ams_data -> Gy);
|
||||
if (ams_data -> dim >= 2)
|
||||
hypre_ParVectorDestroy(ams_data -> Gy);
|
||||
if (ams_data -> dim == 3)
|
||||
hypre_ParVectorDestroy(ams_data -> Gz);
|
||||
}
|
||||
@ -2144,27 +2233,30 @@ HYPRE_Int hypre_AMSSetup(void *solver,
|
||||
(HYPRE_ParCSRMatrix)ams_data -> A_Pix,
|
||||
0, 0);
|
||||
|
||||
if (!hypre_ParCSRMatrixCommPkg(ams_data -> Piy))
|
||||
hypre_MatvecCommPkgCreate(ams_data -> Piy);
|
||||
P_owned_col_starts = hypre_ParCSRMatrixOwnsColStarts(ams_data -> Piy);
|
||||
hypre_BoomerAMGBuildCoarseOperator(ams_data -> Piy,
|
||||
ams_data -> A,
|
||||
ams_data -> Piy,
|
||||
&ams_data -> A_Piy);
|
||||
if (!P_owned_col_starts)
|
||||
if (ams_data -> Piy)
|
||||
{
|
||||
hypre_ParCSRMatrixOwnsRowStarts(ams_data -> A_Piy) = 0;
|
||||
hypre_ParCSRMatrixOwnsColStarts(ams_data -> A_Piy) = 0;
|
||||
if (!hypre_ParCSRMatrixCommPkg(ams_data -> Piy))
|
||||
hypre_MatvecCommPkgCreate(ams_data -> Piy);
|
||||
P_owned_col_starts = hypre_ParCSRMatrixOwnsColStarts(ams_data -> Piy);
|
||||
hypre_BoomerAMGBuildCoarseOperator(ams_data -> Piy,
|
||||
ams_data -> A,
|
||||
ams_data -> Piy,
|
||||
&ams_data -> A_Piy);
|
||||
if (!P_owned_col_starts)
|
||||
{
|
||||
hypre_ParCSRMatrixOwnsRowStarts(ams_data -> A_Piy) = 0;
|
||||
hypre_ParCSRMatrixOwnsColStarts(ams_data -> A_Piy) = 0;
|
||||
}
|
||||
|
||||
/* Make sure that A_Piy has no zero rows (this can happen
|
||||
for some kinds of boundary conditions with contact). */
|
||||
hypre_ParCSRMatrixFixZeroRows(ams_data -> A_Piy);
|
||||
|
||||
HYPRE_BoomerAMGSetup(ams_data -> B_Piy,
|
||||
(HYPRE_ParCSRMatrix)ams_data -> A_Piy,
|
||||
0, 0);
|
||||
}
|
||||
|
||||
/* Make sure that A_Piy has no zero rows (this can happen
|
||||
for some kinds of boundary conditions with contact). */
|
||||
hypre_ParCSRMatrixFixZeroRows(ams_data -> A_Piy);
|
||||
|
||||
HYPRE_BoomerAMGSetup(ams_data -> B_Piy,
|
||||
(HYPRE_ParCSRMatrix)ams_data -> A_Piy,
|
||||
0, 0);
|
||||
|
||||
if (ams_data -> Piz)
|
||||
{
|
||||
if (!hypre_ParCSRMatrixCommPkg(ams_data -> Piz))
|
||||
@ -2278,7 +2370,8 @@ HYPRE_Int hypre_AMSSetup(void *solver,
|
||||
if (input_info == 1)
|
||||
{
|
||||
hypre_ParVectorDestroy(ams_data -> Gx);
|
||||
hypre_ParVectorDestroy(ams_data -> Gy);
|
||||
if (ams_data -> dim >= 2)
|
||||
hypre_ParVectorDestroy(ams_data -> Gy);
|
||||
if (ams_data -> dim == 3)
|
||||
hypre_ParVectorDestroy(ams_data -> Gz);
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user