Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 16 additions & 33 deletions cscs-checks/libraries/petsc/petsc_helloworld.py
Original file line number Diff line number Diff line change
@@ -1,24 +1,13 @@
import filecmp
import os
import reframe as rfm
import reframe.utility.sanity as sn

from reframe.core.pipeline import RegressionTest


@sn.sanity_function
def sanity_filecmp(output, reference):
return sn.assert_true(
filecmp.cmp(output, reference, shallow=False),
msg="files are not the same: `%s', `%s'" % (output, reference)
)


class PetscPoisson2DCheck(RegressionTest):
def __init__(self, variant, **kwargs):
super().__init__('petsc_2dpoisson_%s' % variant,
os.path.dirname(__file__), **kwargs)
@rfm.parameterized_test(['dynamic'], ['static'])
class PetscPoisson2DCheck(rfm.RegressionTest):
def __init__(self, variant):
super().__init__()
self.descr = ('Compile/run PETSc 2D Poisson example with cray-petsc '
'(%s linking case)') % variant
'(%s linking)') % variant
self.valid_systems = ['daint:gpu', 'daint:mc',
'dom:gpu', 'dom:mc']
self.valid_prog_environs = ['PrgEnv-cray', 'PrgEnv-gnu',
Expand All @@ -27,22 +16,16 @@ def __init__(self, variant, **kwargs):
self.modules = ['cray-petsc']
self.num_tasks = 16
self.num_tasks_per_node = 8
self.dynamic = True if variant == 'dynamic' else False
self.executable_opts = ['-da_grid_x 4', '-da_grid_y 4', '-mat_view',
'> petsc_poisson2d.out']
self.build_system = 'SingleSource'
if variant == 'dynamic':
self.build_system.cflags = ['-dynamic']

self.maintainers = ['WS', 'AJ']
self.tags = {'production'}
self.executable_opts = ['-da_grid_x 4', '-da_grid_y 4', '-ksp_monitor']

self.sanity_patterns = sanity_filecmp(
'petsc_poisson2d.out', 'petsc_poisson2d.ref')
# Check the final residual norm for convergence
norm = sn.extractsingle(r'\s+\d+\s+KSP Residual norm\s+(?P<norm>\S+)',
self.stdout, 'norm', float, -1)
self.sanity_patterns = sn.assert_lt(norm, 1.0e-5)

def compile(self):
if self.dynamic:
self.current_environ.cxxflags = '-dynamic'
super().compile()


def _get_checks(**kwargs):
return [PetscPoisson2DCheck('dynamic', **kwargs),
PetscPoisson2DCheck('static', **kwargs)]
self.tags = {'production'}
self.maintainers = ['WS', 'AJ', 'TM']
36 changes: 0 additions & 36 deletions cscs-checks/libraries/petsc/src/petsc_poisson2d.ref

This file was deleted.

116 changes: 58 additions & 58 deletions cscs-checks/libraries/petsc/src/poisson2d.c
Original file line number Diff line number Diff line change
Expand Up @@ -11,11 +11,13 @@
Contributed by Michael Boghosian <boghmic@iit.edu>, 2008,
based on petsc/src/ksp/ksp/examples/tutorials/ex29.c and ex32.c

Compare to ex66.c

Example of Usage:
./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 3 -ksp_monitor -ksp_view -dm_view draw -draw_pause -1
./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type ilu -mg_levels_0_pc_factor_levels 1 -ksp_monitor -ksp_view
./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor
mpiexec -n 4 ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 10 -ksp_monitor -ksp_view -log_summary
mpiexec -n 4 ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 10 -ksp_monitor -ksp_view -log_view
*/

static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";
Expand All @@ -28,37 +30,28 @@ static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";

extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,void*);
extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
extern PetscErrorCode ComputeTrueSolution(DM, Vec);
extern PetscErrorCode VecView_VTK(Vec, const char [], const char []);

typedef enum {DIRICHLET, NEUMANN} BCType;

typedef struct {
PetscScalar uu, tt;
BCType bcType;
} UserContext;

#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **argv)
{
KSP ksp;
DM da;
UserContext user;
PetscInt bc;
PetscErrorCode ierr;

PetscInitialize(&argc,&argv,(char*)0,help);
ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,(DM)da);
ierr = DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,11,11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);CHKERRQ(ierr);
ierr = DMSetFromOptions(da);CHKERRQ(ierr);
ierr = DMSetUp(da);CHKERRQ(ierr);
ierr = KSPSetDM(ksp,(DM)da);CHKERRQ(ierr);
ierr = DMSetApplicationContext(da,&user);CHKERRQ(ierr);

user.uu = 1.0;
user.tt = 1.0;
bc = (PetscInt)NEUMANN; /* Use Neumann Boundary Conditions */
user.bcType = (BCType)bc;


ierr = KSPSetComputeRHS(ksp,ComputeRHS,&user);CHKERRQ(ierr);
ierr = KSPSetComputeOperators(ksp,ComputeJacobian,&user);CHKERRQ(ierr);
Expand All @@ -67,12 +60,10 @@ int main(int argc,char **argv)

ierr = DMDestroy(&da);CHKERRQ(ierr);
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = PetscFinalize();CHKERRQ(ierr);
return 0;
ierr = PetscFinalize();
return ierr;
}

#undef __FUNCT__
#define __FUNCT__ "ComputeRHS"
PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
{
UserContext *user = (UserContext*)ctx;
Expand All @@ -81,6 +72,7 @@ PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
PetscScalar Hx,Hy,pi,uu,tt;
PetscScalar **array;
DM da;
MatNullSpace nullspace;

PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
Expand All @@ -91,7 +83,6 @@ PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
Hy = 1.0/(PetscReal)(N);

ierr = DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);CHKERRQ(ierr); /* Fine grid */
/* printf(" M N: %d %d; xm ym: %d %d; xs ys: %d %d\n",M,N,xm,ym,xs,ys); */
ierr = DMDAVecGetArray(da, b, &array);CHKERRQ(ierr);
for (j=ys; j<ys+ym; j++) {
for (i=xs; i<xs+xm; i++) {
Expand All @@ -104,26 +95,20 @@ PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)

/* force right hand side to be consistent for singular matrix */
/* note this is really a hack, normally the model would provide you with a consistent right handside */
if (user->bcType == NEUMANN) {
MatNullSpace nullspace;

ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
}
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
ierr = MatNullSpaceRemove(nullspace,b);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "ComputeJacobian"
PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
{
UserContext *user = (UserContext*)ctx;
PetscErrorCode ierr;
PetscInt i, j, M, N, xm, ym, xs, ys, num, numi, numj;
PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
MatStencil row, col[5];
DM da;
MatNullSpace nullspace;

PetscFunctionBeginUser;
ierr = KSPGetDM(ksp,&da);CHKERRQ(ierr);
Expand All @@ -138,29 +123,26 @@ PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
row.i = i; row.j = j;

if (i==0 || j==0 || i==M-1 || j==N-1) {
if (user->bcType == DIRICHLET) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
else if (user->bcType == NEUMANN) {
num=0; numi=0; numj=0;
if (j!=0) {
v[num] = -HxdHy; col[num].i = i; col[num].j = j-1;
num++; numj++;
}
if (i!=0) {
v[num] = -HydHx; col[num].i = i-1; col[num].j = j;
num++; numi++;
}
if (i!=M-1) {
v[num] = -HydHx; col[num].i = i+1; col[num].j = j;
num++; numi++;
}
if (j!=N-1) {
v[num] = -HxdHy; col[num].i = i; col[num].j = j+1;
num++; numj++;
}
v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i; col[num].j = j;
num++;
ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
num=0; numi=0; numj=0;
if (j!=0) {
v[num] = -HxdHy; col[num].i = i; col[num].j = j-1;
num++; numj++;
}
if (i!=0) {
v[num] = -HydHx; col[num].i = i-1; col[num].j = j;
num++; numi++;
}
if (i!=M-1) {
v[num] = -HydHx; col[num].i = i+1; col[num].j = j;
num++; numi++;
}
if (j!=N-1) {
v[num] = -HxdHy; col[num].i = i; col[num].j = j+1;
num++; numj++;
}
v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i; col[num].j = j;
num++;
ierr = MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);CHKERRQ(ierr);
} else {
v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
Expand All @@ -173,13 +155,31 @@ PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
}
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (user->bcType == NEUMANN) {
MatNullSpace nullspace;

ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
ierr = MatSetNullSpace(jac,nullspace);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
}
ierr = MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);CHKERRQ(ierr);
ierr = MatSetNullSpace(J,nullspace);CHKERRQ(ierr);
ierr = MatNullSpaceDestroy(&nullspace);CHKERRQ(ierr);
PetscFunctionReturn(0);
}



/*TEST

build:
requires: !complex !single

test:
args: -pc_type mg -pc_mg_type full -ksp_type cg -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type svd -ksp_view

test:
suffix: 2
nsize: 4
args: -pc_type mg -pc_mg_type full -ksp_type cg -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd -ksp_view

test:
suffix: 3
nsize: 2
args: -pc_type mg -pc_mg_type full -ksp_monitor_short -da_refine 5 -mg_coarse_ksp_type cg -mg_coarse_ksp_converged_reason -mg_coarse_ksp_rtol 1e-2 -mg_coarse_ksp_max_it 5 -mg_coarse_pc_type none -pc_mg_levels 2 -ksp_type pipefgmres -ksp_pipefgmres_shift 1.5

TEST*/