CFD/hypre
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Extended ADS to support (arbitrary) high-order Raviart-Thomas discretizations.

Browse Source 
Modified the high-order AMS interface.
This commit is contained in:
kolev1 2011-11-03 23:38:35 +00:00
parent 71b4a5fc1e
commit b9b84bcd46
7 changed files with 259 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
parcsr_ls/HYPRE_ads.c
View File

@ -103,6 +103,32 @@ HYPRE_Int HYPRE_ADSSetCoordinateVectors(HYPRE_Solver solver,
(hypre_ParVector *) z);
}
/*--------------------------------------------------------------------------
* HYPRE_ADSSetInterpolations
*--------------------------------------------------------------------------*/
HYPRE_Int HYPRE_ADSSetInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix RT_Pi,
HYPRE_ParCSRMatrix RT_Pix,
HYPRE_ParCSRMatrix RT_Piy,
HYPRE_ParCSRMatrix RT_Piz,
HYPRE_ParCSRMatrix ND_Pi,
HYPRE_ParCSRMatrix ND_Pix,
HYPRE_ParCSRMatrix ND_Piy,
HYPRE_ParCSRMatrix ND_Piz)
{
return hypre_ADSSetInterpolations((void *) solver,
(hypre_ParCSRMatrix *) RT_Pi,
(hypre_ParCSRMatrix *) RT_Pix,
(hypre_ParCSRMatrix *) RT_Piy,
(hypre_ParCSRMatrix *) RT_Piz,
(hypre_ParCSRMatrix *) ND_Pi,
(hypre_ParCSRMatrix *) ND_Pix,
(hypre_ParCSRMatrix *) ND_Piy,
(hypre_ParCSRMatrix *) ND_Piz);
}
/*--------------------------------------------------------------------------
* HYPRE_ADSSetMaxIter
*--------------------------------------------------------------------------*/

31
parcsr_ls/HYPRE_ams.c
View File

@ -117,29 +117,20 @@ HYPRE_Int HYPRE_AMSSetEdgeConstantVectors(HYPRE_Solver solver,
}
/*--------------------------------------------------------------------------
* HYPRE_AMSSetNedelecInterpolation
* HYPRE_AMSSetInterpolations
*--------------------------------------------------------------------------*/
HYPRE_Int HYPRE_AMSSetNedelecInterpolation(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pi)
HYPRE_Int HYPRE_AMSSetInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pi,
HYPRE_ParCSRMatrix Pix,
HYPRE_ParCSRMatrix Piy,
HYPRE_ParCSRMatrix Piz)
{
return hypre_AMSSetNedelecInterpolation((void *) solver,
(hypre_ParCSRMatrix *) Pi);
}
/*--------------------------------------------------------------------------
* HYPRE_AMSSetNedelecInterpolations
*--------------------------------------------------------------------------*/
HYPRE_Int HYPRE_AMSSetNedelecInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pix,
HYPRE_ParCSRMatrix Piy,
HYPRE_ParCSRMatrix Piz)
{
return hypre_AMSSetNedelecInterpolations((void *) solver,
(hypre_ParCSRMatrix *) Pix,
(hypre_ParCSRMatrix *) Piy,
(hypre_ParCSRMatrix *) Piz);
return hypre_AMSSetInterpolations((void *) solver,
(hypre_ParCSRMatrix *) Pi,
(hypre_ParCSRMatrix *) Pix,
(hypre_ParCSRMatrix *) Piy,
(hypre_ParCSRMatrix *) Piz);
}
/*--------------------------------------------------------------------------

134
parcsr_ls/HYPRE_parcsr_ls.h
View File

@ -1391,15 +1391,19 @@ HYPRE_Int HYPRE_AMSSetEdgeConstantVectors(HYPRE_Solver solver,
HYPRE_ParVector Gz);
/**
* Sets the Nedelec interpolation matrix $\Pi$.
* (Optional) Set the (components of) the Nedelec interpolation matrix
* $\Pi = [ \Pi^x, \Pi^y, \Pi^z ]$.
*
* This is an optional call, which is generally intended to be used only for
* high-order Nedelec discretizations (in the lowest order case, $\Pi$ is
* constructed internally in AMS from the discreet gradient matrix and the
* coordinates of the vertices). By definition, $\Pi$ is the matrix
* representation of the linear operator that interpolates (high-order) vector
* nodal finite elements into the (high-order) Nedelec space. Note that this
* depends on the choice for the basis in the (high-order) nodal space.
* This function is generally intended to be used only for high-order Nedelec
* discretizations (in the lowest order case, $\Pi$ is constructed internally in
* AMS from the discreet gradient matrix and the coordinates of the vertices).
*
* By definition, $\Pi$ is the matrix representation of the linear operator that
* interpolates (high-order) vector nodal finite elements into the (high-order)
* Nedelec space. The component matrices are defined as $\Pi^x \varphi = \Pi
* (\varphi,0,0)$ and similarly for $\Pi^y$ and $\Pi^z$. Note that all these
* operators depend on the choice of the basis and degrees of freedom in the
* high-order spaces.
*
* The column numbering of Pi should be node-based, i.e. the $x$/$y$/$z$
* components of the first node (vertex or high-order dof) should be listed
@ -1407,33 +1411,17 @@ HYPRE_Int HYPRE_AMSSetEdgeConstantVectors(HYPRE_Solver solver,
* (see the documentation of HYPRE\_BoomerAMGSetDofFunc).
*
* If used, this function should be called before HYPRE\_AMSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only AMS cycle types
* based on monolithic $\Pi$ (i.e. cycle\_type < 10) will be available.
* no need to provide the vertex coordinates. Furthermore, only one of the sets
* $\{\Pi\}$ and $\{\Pi^x,\Pi^y,\Pi^z\}$ needs to be specified (though it is OK
* to provide both). If Pix is NULL, then scalar $\Pi$-based AMS cycles,
* i.e. those with cycle\_type > 10, will be unavailable. Similarly, AMS cycles
* based on monolithic $\Pi$ (cycle\_type < 10) require that Pi is not NULL.
**/
HYPRE_Int HYPRE_AMSSetNedelecInterpolation(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pi);
/**
* Set the components of the Nedelec interpolation matrix $\Pi = [ \Pi^x, \Pi^y, \Pi^z ]$.
*
* This is an optional call, which is generally intended to be used only for
* high-order Nedelec discretizations (in the lowest order case, $\Pi$ is
* constructed internally in AMS from the discrete gradient matrix and the
* coordinates of the vertices). By definition, $\Pi^x$ is the matrix
* representation of the linear operator that interpolates (high-order) vector
* nodal finite elements in the form $(\varphi,0,0)$ into the (high-order) Nedelec
* space. In other words, $\Pi^x \varphi = \Pi (\varphi,0,0)$ and similarly for
* $\Pi^y$ and $\Pi^z$. Note that these depends on the choice for the basis in
* the (high-order) nodal space.
*
* If used, this function should be called before HYPRE\_AMSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only AMS cycle types
* based on scalar $\Pi$ (i.e. cycle\_type > 10) will be available.
**/
HYPRE_Int HYPRE_AMSSetNedelecInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pix,
HYPRE_ParCSRMatrix Piy,
HYPRE_ParCSRMatrix Piz);
HYPRE_Int HYPRE_AMSSetInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix Pi,
HYPRE_ParCSRMatrix Pix,
HYPRE_ParCSRMatrix Piy,
HYPRE_ParCSRMatrix Piz);
/**
* (Optional) Sets the matrix $A_\alpha$ corresponding to the Poisson
@ -1629,12 +1617,12 @@ HYPRE_Int HYPRE_AMSConstructDiscreteGradient(HYPRE_ParCSRMatrix A,
/**
* Create an ADS solver object.
**/
HYPRE_Int HYPRE_ADSCreate ( HYPRE_Solver *solver );
HYPRE_Int HYPRE_ADSCreate(HYPRE_Solver *solver);
/**
* Destroy an ADS solver object.
**/
HYPRE_Int HYPRE_ADSDestroy ( HYPRE_Solver solver );
HYPRE_Int HYPRE_ADSDestroy(HYPRE_Solver solver);
/**
* Set up the ADS solver or preconditioner.
@ -1646,7 +1634,7 @@ HYPRE_Int HYPRE_ADSDestroy ( HYPRE_Solver solver );
* @param b Ignored by this function.
* @param x Ignored by this function.
**/
HYPRE_Int HYPRE_ADSSetup ( HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_ParVector b , HYPRE_ParVector x );
HYPRE_Int HYPRE_ADSSetup(HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_ParVector b , HYPRE_ParVector x);
/**
* Solve the system or apply ADS as a preconditioner.
@ -1658,40 +1646,82 @@ HYPRE_Int HYPRE_ADSSetup ( HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_Pa
* @param b [IN] right hand side of the linear system to be solved
* @param x [OUT] approximated solution of the linear system to be solved
**/
HYPRE_Int HYPRE_ADSSolve ( HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_ParVector b , HYPRE_ParVector x );
HYPRE_Int HYPRE_ADSSolve(HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_ParVector b , HYPRE_ParVector x);
/**
* Sets the discrete curl matrix $C$.
* This function should be called before HYPRE\_ADSSetup()!
**/
HYPRE_Int HYPRE_ADSSetDiscreteCurl ( HYPRE_Solver solver , HYPRE_ParCSRMatrix C );
HYPRE_Int HYPRE_ADSSetDiscreteCurl(HYPRE_Solver solver , HYPRE_ParCSRMatrix C);
/**
* Sets the discrete gradient matrix $G$.
* This function should be called before HYPRE\_ADSSetup()!
**/
HYPRE_Int HYPRE_ADSSetDiscreteGradient ( HYPRE_Solver solver , HYPRE_ParCSRMatrix G );
HYPRE_Int HYPRE_ADSSetDiscreteGradient(HYPRE_Solver solver , HYPRE_ParCSRMatrix G);
/**
* Sets the $x$, $y$ and $z$ coordinates of the vertices in the mesh.
* This function should be called before HYPRE\_ADSSetup()!
**/
HYPRE_Int HYPRE_ADSSetCoordinateVectors ( HYPRE_Solver solver , HYPRE_ParVector x , HYPRE_ParVector y ,
HYPRE_ParVector z );
HYPRE_Int HYPRE_ADSSetCoordinateVectors(HYPRE_Solver solver , HYPRE_ParVector x , HYPRE_ParVector y ,
HYPRE_ParVector z);
/**
* (Optional) Set the (components of) the Raviart-Thomas ($\Pi_{RT}$) and the Nedelec
* ($\Pi_{ND}$) interpolation matrices.
*
* This function is generally intended to be used only for high-order $H(div)$
* discretizations (in the lowest order case, these matrices are constructed
* internally in ADS from the discreet gradient and curl matrices and the
* coordinates of the vertices).
*
* By definition, RT\_Pi and ND\_Pi are the matrix representations of the linear
* operators $\Pi_{RT}$ and $\Pi_{ND}$ that interpolate (high-order) vector
* nodal finite elements into the (high-order) Raviart-Thomas and Nedelec
* spaces. The component matrices are defined in both cases as $\Pi^x \varphi =
* \Pi (\varphi,0,0)$ and similarly for $\Pi^y$ and $\Pi^z$. Note that all these
* operators depend on the choice of the basis and degrees of freedom in the
* high-order spaces.
*
* The column numbering of RT\_Pi and ND\_Pi should be node-based, i.e. the
* $x$/$y$/$z$ components of the first node (vertex or high-order dof) should be
* listed first, followed by the $x$/$y$/$z$ components of the second node and
* so on (see the documentation of HYPRE\_BoomerAMGSetDofFunc).
*
* If used, this function should be called before hypre\_ADSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only one of the sets
* $\{\Pi_{RT}\}$ and $\{\Pi_{RT}^x,\Pi_{RT}^y,\Pi_{RT}^z\}$ needs to be
* specified (though it is OK to provide both). If RT\_Pix is NULL, then scalar
* $\Pi$-based ADS cycles, i.e. those with cycle\_type > 10, will be
* unavailable. Similarly, ADS cycles based on monolithic $\Pi$ (cycle\_type <
* 10) require that RT\_Pi is not NULL. The same restrictions hold for the sets
* $\{\Pi_{ND}\}$ and $\{\Pi_{ND}^x,\Pi_{ND}^y,\Pi_{ND}^z\}$ -- only one of them
* needs to be specified, and the availability of each enables different AMS
* cycle type options.
**/
HYPRE_Int HYPRE_ADSSetInterpolations(HYPRE_Solver solver,
HYPRE_ParCSRMatrix RT_Pi,
HYPRE_ParCSRMatrix RT_Pix,
HYPRE_ParCSRMatrix RT_Piy,
HYPRE_ParCSRMatrix RT_Piz,
HYPRE_ParCSRMatrix ND_Pi,
HYPRE_ParCSRMatrix ND_Pix,
HYPRE_ParCSRMatrix ND_Piy,
HYPRE_ParCSRMatrix ND_Piz);
/**
* (Optional) Sets maximum number of iterations, if ADS is used
* as a solver. To use ADS as a preconditioner, set the maximum
* number of iterations to $1$. The default is $20$.
**/
HYPRE_Int HYPRE_ADSSetMaxIter ( HYPRE_Solver solver , HYPRE_Int maxit );
HYPRE_Int HYPRE_ADSSetMaxIter(HYPRE_Solver solver , HYPRE_Int maxit);
/**
* (Optional) Set the convergence tolerance, if ADS is used
* as a solver. When using ADS as a preconditioner, set the tolerance
* to $0.0$. The default is $10^{-6}$.
**/
HYPRE_Int HYPRE_ADSSetTol ( HYPRE_Solver solver , double tol );
HYPRE_Int HYPRE_ADSSetTol(HYPRE_Solver solver , double tol);
/**
* (Optional) Choose which auxiliary-space solver to use. Possible values are:
@ -1715,14 +1745,14 @@ HYPRE_Int HYPRE_ADSSetTol ( HYPRE_Solver solver , double tol );
*
* The default is $1$. See the user's manual for more details.
**/
HYPRE_Int HYPRE_ADSSetCycleType ( HYPRE_Solver solver , HYPRE_Int cycle_type );
HYPRE_Int HYPRE_ADSSetCycleType(HYPRE_Solver solver , HYPRE_Int cycle_type);
/**
* (Optional) Control how much information is printed during the
* solution iterations.
* The default is $1$ (print residual norm at each step).
**/
HYPRE_Int HYPRE_ADSSetPrintLevel ( HYPRE_Solver solver , HYPRE_Int print_level );
HYPRE_Int HYPRE_ADSSetPrintLevel(HYPRE_Solver solver , HYPRE_Int print_level);
/**
* (Optional) Sets relaxation parameters for $A$.
@ -1739,35 +1769,35 @@ HYPRE_Int HYPRE_ADSSetPrintLevel ( HYPRE_Solver solver , HYPRE_Int print_level )
* \hline
* \end{tabular}
**/
HYPRE_Int HYPRE_ADSSetSmoothingOptions ( HYPRE_Solver solver , HYPRE_Int relax_type , HYPRE_Int relax_times , double relax_weight , double omega );
HYPRE_Int HYPRE_ADSSetSmoothingOptions(HYPRE_Solver solver , HYPRE_Int relax_type , HYPRE_Int relax_times , double relax_weight , double omega);
/**
* (Optional) Sets parameters for Chebyshev relaxation.
* The defaults are $2$, $0.3$.
**/
HYPRE_Int HYPRE_ADSSetChebySmoothingOptions ( HYPRE_Solver solver , HYPRE_Int cheby_order , HYPRE_Int cheby_fraction );
HYPRE_Int HYPRE_ADSSetChebySmoothingOptions(HYPRE_Solver solver , HYPRE_Int cheby_order , HYPRE_Int cheby_fraction);
/**
* (Optional) Sets AMS parameters for $B_C$.
* The defaults are $11$, $10$, $1$, $3$, $0.25$, $0$, $0$. See the user's manual for more details.
**/
HYPRE_Int HYPRE_ADSSetAMSOptions ( HYPRE_Solver solver , HYPRE_Int cycle_type , HYPRE_Int coarsen_type , HYPRE_Int agg_levels , HYPRE_Int relax_type , double strength_threshold , HYPRE_Int interp_type , HYPRE_Int Pmax );
HYPRE_Int HYPRE_ADSSetAMSOptions(HYPRE_Solver solver , HYPRE_Int cycle_type , HYPRE_Int coarsen_type , HYPRE_Int agg_levels , HYPRE_Int relax_type , double strength_threshold , HYPRE_Int interp_type , HYPRE_Int Pmax);
/**
* (Optional) Sets AMG parameters for $B_\Pi$.
* The defaults are $10$, $1$, $3$, $0.25$, $0$, $0$. See the user's manual for more details.
**/
HYPRE_Int HYPRE_ADSSetAMGOptions ( HYPRE_Solver solver , HYPRE_Int coarsen_type , HYPRE_Int agg_levels , HYPRE_Int relax_type , double strength_threshold , HYPRE_Int interp_type , HYPRE_Int Pmax );
HYPRE_Int HYPRE_ADSSetAMGOptions(HYPRE_Solver solver , HYPRE_Int coarsen_type , HYPRE_Int agg_levels , HYPRE_Int relax_type , double strength_threshold , HYPRE_Int interp_type , HYPRE_Int Pmax);
/**
* Returns the number of iterations taken.
**/
HYPRE_Int HYPRE_ADSGetNumIterations ( HYPRE_Solver solver , HYPRE_Int *num_iterations );
HYPRE_Int HYPRE_ADSGetNumIterations(HYPRE_Solver solver , HYPRE_Int *num_iterations);
/**
* Returns the norm of the final relative residual.
**/
HYPRE_Int HYPRE_ADSGetFinalRelativeResidualNorm ( HYPRE_Solver solver , double *rel_resid_norm );
HYPRE_Int HYPRE_ADSGetFinalRelativeResidualNorm(HYPRE_Solver solver , double *rel_resid_norm);
/*@}*/

4
parcsr_ls/_hypre_parcsr_ls.h
View File

@ -512,8 +512,7 @@ HYPRE_Int hypre_AMSSetDimension ( void *solver , HYPRE_Int dim );
HYPRE_Int hypre_AMSSetDiscreteGradient ( void *solver , hypre_ParCSRMatrix *G );
HYPRE_Int hypre_AMSSetCoordinateVectors ( void *solver , hypre_ParVector *x , hypre_ParVector *y , hypre_ParVector *z );
HYPRE_Int hypre_AMSSetEdgeConstantVectors ( void *solver , hypre_ParVector *Gx , hypre_ParVector *Gy , hypre_ParVector *Gz );
HYPRE_Int hypre_AMSSetNedelecInterpolation ( void *solver , hypre_ParCSRMatrix *Pi);
HYPRE_Int hypre_AMSSetNedelecInterpolations ( void *solver , hypre_ParCSRMatrix *Pix , hypre_ParCSRMatrix *Piy , hypre_ParCSRMatrix *Piz);
HYPRE_Int hypre_AMSSetInterpolations ( void *solver , hypre_ParCSRMatrix *Pi , hypre_ParCSRMatrix *Pix , hypre_ParCSRMatrix *Piy , hypre_ParCSRMatrix *Piz);
HYPRE_Int hypre_AMSSetAlphaPoissonMatrix ( void *solver , hypre_ParCSRMatrix *A_Pi );
HYPRE_Int hypre_AMSSetBetaPoissonMatrix ( void *solver , hypre_ParCSRMatrix *A_G );
HYPRE_Int hypre_AMSSetInteriorNodes ( void *solver , hypre_ParVector *interior_nodes );
@ -584,6 +583,7 @@ HYPRE_Int HYPRE_ADSSolve ( HYPRE_Solver solver , HYPRE_ParCSRMatrix A , HYPRE_Pa
HYPRE_Int HYPRE_ADSSetDiscreteCurl ( HYPRE_Solver solver , HYPRE_ParCSRMatrix C );
HYPRE_Int HYPRE_ADSSetDiscreteGradient ( HYPRE_Solver solver , HYPRE_ParCSRMatrix G );
HYPRE_Int HYPRE_ADSSetCoordinateVectors ( HYPRE_Solver solver , HYPRE_ParVector x , HYPRE_ParVector y , HYPRE_ParVector z );
HYPRE_Int hypre_ADSSetInterpolations ( void *solver , hypre_ParCSRMatrix *RT_Pi , hypre_ParCSRMatrix *RT_Pix , hypre_ParCSRMatrix *RT_Piy , hypre_ParCSRMatrix *RT_Piz , hypre_ParCSRMatrix *ND_Pi , hypre_ParCSRMatrix *ND_Pix , hypre_ParCSRMatrix *ND_Piy , hypre_ParCSRMatrix *ND_Piz );
HYPRE_Int HYPRE_ADSSetMaxIter ( HYPRE_Solver solver , HYPRE_Int maxit );
HYPRE_Int HYPRE_ADSSetTol ( HYPRE_Solver solver , double tol );
HYPRE_Int HYPRE_ADSSetCycleType ( HYPRE_Solver solver , HYPRE_Int cycle_type );

142
parcsr_ls/ads.c
View File

@ -92,6 +92,12 @@ void * hypre_ADSCreate()
ads_data -> A_max_eig_est = 0;
ads_data -> A_min_eig_est = 0;
ads_data -> owns_Pi = 1;
ads_data -> ND_Pi = NULL;
ads_data -> ND_Pix = NULL;
ads_data -> ND_Piy = NULL;
ads_data -> ND_Piz = NULL;
return (void *) ads_data;
}
@ -111,26 +117,26 @@ HYPRE_Int hypre_ADSDestroy(void *solver)
if (ads_data -> B_C)
HYPRE_AMSDestroy(ads_data -> B_C);
if (ads_data -> Pi)
if (ads_data -> owns_Pi && ads_data -> Pi)
hypre_ParCSRMatrixDestroy(ads_data -> Pi);
if (ads_data -> A_Pi)
hypre_ParCSRMatrixDestroy(ads_data -> A_Pi);
if (ads_data -> B_Pi)
HYPRE_BoomerAMGDestroy(ads_data -> B_Pi);
if (ads_data -> Pix)
if (ads_data -> owns_Pi && ads_data -> Pix)
hypre_ParCSRMatrixDestroy(ads_data -> Pix);
if (ads_data -> A_Pix)
hypre_ParCSRMatrixDestroy(ads_data -> A_Pix);
if (ads_data -> B_Pix)
HYPRE_BoomerAMGDestroy(ads_data -> B_Pix);
if (ads_data -> Piy)
if (ads_data -> owns_Pi && ads_data -> Piy)
hypre_ParCSRMatrixDestroy(ads_data -> Piy);
if (ads_data -> A_Piy)
hypre_ParCSRMatrixDestroy(ads_data -> A_Piy);
if (ads_data -> B_Piy)
HYPRE_BoomerAMGDestroy(ads_data -> B_Piy);
if (ads_data -> Piz)
if (ads_data -> owns_Pi && ads_data -> Piz)
hypre_ParCSRMatrixDestroy(ads_data -> Piz);
if (ads_data -> A_Piz)
hypre_ParCSRMatrixDestroy(ads_data -> A_Piz);
@ -210,6 +216,63 @@ HYPRE_Int hypre_ADSSetCoordinateVectors(void *solver,
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_ADSSetInterpolations
*
* Set the (components of) the Raviart-Thomas (RT_Pi) and the Nedelec (ND_Pi)
* interpolation matrices.
*
* This function is generally intended to be used only for high-order H(div)
* discretizations (in the lowest order case, these matrices are constructed
* internally in ADS from the discreet gradient and curl matrices and the
* coordinates of the vertices).
*
* By definition, RT_Pi and ND_Pi are the matrix representations of the linear
* operators that interpolate (high-order) vector nodal finite elements into the
* (high-order) Raviart-Thomas and Nedelec spaces. The component matrices are
* defined in both cases as Pix phi = Pi (phi,0,0) and similarly for Piy and
* Piz. Note that all these operators depend on the choice of the basis and
* degrees of freedom in the high-order spaces.
*
* The column numbering of RT_Pi and ND_pi should be node-based, i.e. the x/y/z
* components of the first node (vertex or high-order dof) should be listed
* first, followed by the x/y/z components of the second node and so on (see the
* documentation of HYPRE_BoomerAMGSetDofFunc).
*
* If used, this function should be called before hypre_ADSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only one of the sets
* {RT_Pi} and {RT_Pix,RT_Piy,RT_Piz} needs to be specified (though it is OK to
* provide both). If RT_Pix is NULL, then scalar Pi-based ADS cycles, i.e.
* those with cycle_type > 10, will be unavailable. Similarly, ADS cycles based
* on monolithic Pi (cycle_type < 10) require that RT_Pi is not NULL. The same
* restrictions hold for the sets {ND_Pi} and {ND_Pix,ND_Piy,ND_Piz} -- only one
* of them needs to be specified, and the availability of each enables different
* AMS cycle type options.
*--------------------------------------------------------------------------*/
HYPRE_Int hypre_ADSSetInterpolations(void *solver,
hypre_ParCSRMatrix *RT_Pi,
hypre_ParCSRMatrix *RT_Pix,
hypre_ParCSRMatrix *RT_Piy,
hypre_ParCSRMatrix *RT_Piz,
hypre_ParCSRMatrix *ND_Pi,
hypre_ParCSRMatrix *ND_Pix,
hypre_ParCSRMatrix *ND_Piy,
hypre_ParCSRMatrix *ND_Piz)
{
hypre_ADSData *ads_data = solver;
ads_data -> Pi = RT_Pi;
ads_data -> Pix = RT_Pix;
ads_data -> Piy = RT_Piy;
ads_data -> Piz = RT_Piz;
ads_data -> ND_Pi = ND_Pi;
ads_data -> ND_Pix = ND_Pix;
ads_data -> ND_Piy = ND_Piy;
ads_data -> ND_Piz = ND_Piz;
ads_data -> owns_Pi = 0;
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_ADSSetMaxIter
*
@ -349,10 +412,6 @@ HYPRE_Int hypre_ADSSetAMSOptions(void *solver,
ads_data -> B_C_interp_type = B_C_interp_type;
ads_data -> B_C_Pmax = B_C_Pmax;
if (B_C_cycle_type < 10)
hypre_error_w_msg(HYPRE_ERROR_GENERIC,
"ADS requires AMS solver type greater than 10!");
return hypre_error_flag;
}
@ -843,10 +902,29 @@ HYPRE_Int hypre_ADSSetup(void *solver,
HYPRE_AMSSetCycleType(ads_data -> B_C, ads_data -> B_C_cycle_type);
HYPRE_AMSSetDiscreteGradient(ads_data -> B_C,
(HYPRE_ParCSRMatrix) ads_data -> G);
HYPRE_AMSSetCoordinateVectors(ads_data -> B_C,
(HYPRE_ParVector) ads_data -> x,
(HYPRE_ParVector) ads_data -> y,
(HYPRE_ParVector) ads_data -> z);
if (ads_data -> ND_Pi == NULL && ads_data -> ND_Pix == NULL)
{
if (ads_data -> B_C_cycle_type < 10)
hypre_error_w_msg(HYPRE_ERROR_GENERIC,
"Unsupported AMS cycle type in ADS!");
HYPRE_AMSSetCoordinateVectors(ads_data -> B_C,
(HYPRE_ParVector) ads_data -> x,
(HYPRE_ParVector) ads_data -> y,
(HYPRE_ParVector) ads_data -> z);
}
else
{
if ((ads_data -> B_C_cycle_type < 10 && ads_data -> ND_Pi == NULL) ||
(ads_data -> B_C_cycle_type > 10 && ads_data -> ND_Pix == NULL))
hypre_error_w_msg(HYPRE_ERROR_GENERIC,
"Unsupported AMS cycle type in ADS!");
HYPRE_AMSSetInterpolations(ads_data -> B_C,
(HYPRE_ParCSRMatrix) ads_data -> ND_Pi,
(HYPRE_ParCSRMatrix) ads_data -> ND_Pix,
(HYPRE_ParCSRMatrix) ads_data -> ND_Piy,
(HYPRE_ParCSRMatrix) ads_data -> ND_Piz);
}
/* beta=0 in the subspace */
HYPRE_AMSSetBetaPoissonMatrix(ads_data -> B_C, NULL);
@ -890,9 +968,26 @@ HYPRE_Int hypre_ADSSetup(void *solver,
}
ams_data = (hypre_AMSData *) ads_data -> B_C;
if (ads_data -> cycle_type > 10)
/* Construct Pi{x,y,z} instead of Pi = [Pix,Piy,Piz] */
hypre_ADSComputePixyz(ads_data -> A,
if (ads_data -> Pi == NULL && ads_data -> Pix == NULL)
{
if (ads_data -> cycle_type > 10)
/* Construct Pi{x,y,z} instead of Pi = [Pix,Piy,Piz] */
hypre_ADSComputePixyz(ads_data -> A,
ads_data -> C,
ads_data -> G,
ads_data -> x,
ads_data -> y,
ads_data -> z,
ams_data -> Pix,
ams_data -> Piy,
ams_data -> Piz,
&ads_data -> Pix,
&ads_data -> Piy,
&ads_data -> Piz);
else
/* Construct the Pi interpolation matrix */
hypre_ADSComputePi(ads_data -> A,
ads_data -> C,
ads_data -> G,
ads_data -> x,
@ -901,21 +996,8 @@ HYPRE_Int hypre_ADSSetup(void *solver,
ams_data -> Pix,
ams_data -> Piy,
ams_data -> Piz,
&ads_data -> Pix,
&ads_data -> Piy,
&ads_data -> Piz);
else
/* Construct the Pi interpolation matrix */
hypre_ADSComputePi(ads_data -> A,
ads_data -> C,
ads_data -> G,
ads_data -> x,
ads_data -> y,
ads_data -> z,
ams_data -> Pix,
ams_data -> Piy,
ams_data -> Piz,
&ads_data -> Pi);
&ads_data -> Pi);
}
if (ads_data -> cycle_type > 10)
/* Create the AMG solvers on the range of Pi{x,y,z}^T */

5
parcsr_ls/ads.h
View File

@ -46,6 +46,11 @@ typedef struct
/* AMG solvers for A_Pi{x,y,z} */
HYPRE_Solver B_Pix, B_Piy, B_Piz;
/* Does the solver own the RT/ND interpolations matrices? */
HYPRE_Int owns_Pi;
/* The (high-order) edge interpolation matrix and its components */
hypre_ParCSRMatrix *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz;
/* Discrete gradient matrix (vertex-to-edge) */
hypre_ParCSRMatrix *G;
/* Coordinates of the vertices */

63
parcsr_ls/ams.c
View File

@ -1003,17 +1003,19 @@ HYPRE_Int hypre_AMSSetEdgeConstantVectors(void *solver,
}
/*--------------------------------------------------------------------------
* hypre_AMSSetNedelecInterpolation
* hypre_AMSSetInterpolations
*
* Set the Nedelec interpolation matrix Pi.
* Set the (components of) the Nedelec interpolation matrix Pi=[Pix,Piy,Piz].
*
* This is an optional call, which is generally intended to be used only for
* high-order Nedelec discretizations (in the lowest order case, Pi is
* constructed internally in AMS from the discreet gradient matrix and the
* coordinates of the vertices). By definition, Pi is the matrix representation
* of the linear operator that interpolates (high-order) vector nodal finite
* elements into the (high-order) Nedelec space. Note that this depends on the
* choice for the basis in the (high-order) nodal space.
* This function is generally intended to be used only for high-order Nedelec
* discretizations (in the lowest order case, Pi is constructed internally in
* AMS from the discreet gradient matrix and the coordinates of the vertices).
*
* By definition, Pi is the matrix representation of the linear operator that
* interpolates (high-order) vector nodal finite elements into the (high-order)
* Nedelec space. The component matrices are defined as Pix phi = Pi (phi,0,0)
* and similarly for Piy and Piz. Note that all these operators depend on the
* choice of the basis and degrees of freedom in the high-order spaces.
*
* The column numbering of Pi should be node-based, i.e. the x/y/z components of
* the first node (vertex or high-order dof) should be listed first, followed by
@ -1021,44 +1023,21 @@ HYPRE_Int hypre_AMSSetEdgeConstantVectors(void *solver,
* HYPRE_BoomerAMGSetDofFunc).
*
* If used, this function should be called before hypre_AMSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only AMS cycle types
* based on monolithic Pi (i.e. cycle_type < 10) will be available.
* no need to provide the vertex coordinates. Furthermore, only one of the sets
* {Pi} and {Pix,Piy,Piz} needs to be specified (though it is OK to provide
* both). If Pix is NULL, then scalar Pi-based AMS cycles, i.e. those with
* cycle_type > 10, will be unavailable. Similarly, AMS cycles based on
* monolithic Pi (cycle_type < 10) require that Pi is not NULL.
*--------------------------------------------------------------------------*/
HYPRE_Int hypre_AMSSetNedelecInterpolation(void *solver,
hypre_ParCSRMatrix *Pi)
HYPRE_Int hypre_AMSSetInterpolations(void *solver,
hypre_ParCSRMatrix *Pi,
hypre_ParCSRMatrix *Pix,
hypre_ParCSRMatrix *Piy,
hypre_ParCSRMatrix *Piz)
{
hypre_AMSData *ams_data = solver;
ams_data -> Pi = Pi;
ams_data -> owns_Pi = 0;
return hypre_error_flag;
}
/*--------------------------------------------------------------------------
* hypre_AMSSetNedelecInterpolations
*
* Set the components of the Nedelec interpolation matrix Pi = [Pix, Piy, Piz].
*
* This is an optional call, which is generally intended to be used only for
* high-order Nedelec discretizations (in the lowest order case, Pi is
* constructed internally in AMS from the discrete gradient matrix and the
* coordinates of the vertices). By definition, Pix is the matrix representation
* of the linear operator that interpolates (high-order) vector nodal finite
* elements in the form (phi,0,0) into the (high-order) Nedelec space. In other
* words, Pix phi = Pi (phi,0,0) and similarly for Piy and Piz. Note that these
* depends on the choice for the basis in the (high-order) nodal space.
*
* If used, this function should be called before hypre_AMSSetup() and there is
* no need to provide the vertex coordinates. Furthermore, only AMS cycle types
* based on scalar Pi (i.e. cycle_type > 10) will be available.
*--------------------------------------------------------------------------*/
HYPRE_Int hypre_AMSSetNedelecInterpolations(void *solver,
hypre_ParCSRMatrix *Pix,
hypre_ParCSRMatrix *Piy,
hypre_ParCSRMatrix *Piz)
{
hypre_AMSData *ams_data = solver;
ams_data -> Pix = Pix;
ams_data -> Piy = Piy;
ams_data -> Piz = Piz;