Skip to content

Commit 71b1bc8

Browse files
rezgarshakerirezgarshakeri
committed
Added setup-solver.c
1 parent 4af980f commit 71b1bc8

File tree

9 files changed

+136
-87
lines changed

9 files changed

+136
-87
lines changed

‎examples/Hdiv-mixed/conv_plot.py‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,9 @@
1616
# software, applications, hardware, advanced system engineering and early
1717
# testbed platforms, in support of the nation's exascale computing imperative.
1818

19+
# After ./conv_test.sh you can plot using
20+
# python conv_plot.py -f conv_test_result.csv
21+
1922
import pandas as pd
2023
import argparse
2124
from pylab import *

‎examples/Hdiv-mixed/conv_test.sh‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,9 @@
1616
# software, applications, hardware, advanced system engineering and early
1717
# testbed platforms, in support of the nation's exascale computing imperative.
1818

19-
# After make the problem, you can run convergence test by: ./conv_test.sh -d 2 (or -d 3)
19+
# After make the problem, you can run convergence test by:
20+
#./conv_test.sh -d 2 (or -d 3)
21+
2022
# Reading arguments with getopts options
2123
while getopts d: flag
2224
do

‎examples/Hdiv-mixed/include/setup-libceed.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@ PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P,
1818
CeedElemRestriction *elem_restr_oriented, CeedElemRestriction *elem_restr);
1919
// Set up libCEED for a given degree
2020
PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
21-
ProblemData *problem_data, PetscInt U_g_size,
22-
PetscInt U_loc_size, CeedData ceed_data,
21+
ProblemData *problem_data,
22+
PetscInt rhs_loc_size, CeedData ceed_data,
2323
CeedVector rhs_ceed, CeedVector *target);
2424
#endif // setuplibceed_h

‎examples/Hdiv-mixed/include/setup-solver.h‎

Lines changed: 18 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,18 @@
1+
#ifndef setupsolver_h
2+
#define setupsolver_h
3+
4+
#include <petsc.h>
5+
#include <petscdmplex.h>
6+
#include <petscsys.h>
7+
#include <ceed.h>
8+
#include "../include/structs.h"
9+
#include "../include/matops.h"
10+
11+
// ---------------------------------------------------------------------------
12+
// Setup Solver
13+
// ---------------------------------------------------------------------------
14+
PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed,
15+
CeedData ceed_data, User user, VecType vec_type,
16+
CeedVector target, Vec rhs, Vec *U_g);
17+
18+
#endif // setupsolver_h

‎examples/Hdiv-mixed/main.c‎

Lines changed: 18 additions & 79 deletions
Original file line numberDiff line numberDiff line change
@@ -110,104 +110,47 @@ int main(int argc, char **argv) {
110110
}
111111
PetscCall( DMSetVecType(dm, vec_type) );
112112
}
113-
// ---------------------------------------------------------------------------
114-
// Create global, local solution, local rhs vector
115-
// ---------------------------------------------------------------------------
116-
Vec U_g, U_loc;
117-
PetscInt U_l_size, U_g_size, U_loc_size;
118-
// Create global and local solution vectors
119-
PetscCall( DMCreateGlobalVector(dm, &U_g) );
120-
PetscCall( VecGetSize(U_g, &U_g_size) );
121-
// Local size for matShell
122-
PetscCall( VecGetLocalSize(U_g, &U_l_size) );
123-
// Create local unknown vector U_loc
124-
PetscCall( DMCreateLocalVector(dm, &U_loc) );
125-
// Local size for libCEED
126-
PetscCall( VecGetSize(U_loc, &U_loc_size) );
127113

128-
// Get RHS vector
114+
// ---------------------------------------------------------------------------
115+
// Create local RHS vector
116+
// ---------------------------------------------------------------------------
129117
Vec rhs_loc;
118+
PetscInt rhs_loc_size;
130119
PetscScalar *r;
131120
CeedVector rhs_ceed, target;
132121
PetscMemType mem_type;
133-
PetscCall( VecDuplicate(U_loc, &rhs_loc) );
122+
PetscCall( DMCreateLocalVector(dm, &rhs_loc) );
123+
// Local size for libCEED
124+
PetscCall( VecGetSize(rhs_loc, &rhs_loc_size) );
134125
PetscCall( VecZeroEntries(rhs_loc) );
135126
PetscCall( VecGetArrayAndMemType(rhs_loc, &r, &mem_type) );
136-
CeedVectorCreate(ceed, U_l_size, &rhs_ceed);
127+
CeedVectorCreate(ceed, rhs_loc_size, &rhs_ceed);
137128
CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r);
138129

139130
// ---------------------------------------------------------------------------
140-
// Setup libCEED
131+
// Setup libCEED - Compute local rhs and true solution (target)
141132
// ---------------------------------------------------------------------------
142133
// -- Set up libCEED objects
143-
PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data, U_g_size,
144-
U_loc_size, ceed_data, rhs_ceed, &target) );
134+
PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
135+
rhs_loc_size, ceed_data, rhs_ceed, &target) );
145136
//CeedVectorView(rhs_ceed, "%12.8f", stdout);
137+
146138
// ---------------------------------------------------------------------------
147-
// Gather RHS
139+
// Create global rhs
148140
// ---------------------------------------------------------------------------
149141
Vec rhs;
150142
CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL);
151143
PetscCall( VecRestoreArrayAndMemType(rhs_loc, &r) );
152-
PetscCall( VecDuplicate(U_g, &rhs) );
144+
PetscCall( DMCreateGlobalVector(dm, &rhs) );
153145
PetscCall( VecZeroEntries(rhs) );
154146
PetscCall( DMLocalToGlobal(dm, rhs_loc, ADD_VALUES, rhs) );
155-
// ---------------------------------------------------------------------------
156-
// Setup Mat, KSP
157-
// ---------------------------------------------------------------------------
158-
user->comm = comm;
159-
user->dm = dm;
160-
user->X_loc = U_loc;
161-
PetscCall( VecDuplicate(U_loc, &user->Y_loc) );
162-
user->x_ceed = ceed_data->x_ceed;
163-
user->y_ceed = ceed_data->y_ceed;
164-
user->op_apply = ceed_data->op_residual;
165-
user->op_error = ceed_data->op_error;
166-
user->elem_restr_u = ceed_data->elem_restr_u;
167-
user->ceed = ceed;
168-
// Operator
169-
Mat mat;
170-
PetscCall( MatCreateShell(comm, U_l_size, U_l_size, U_g_size, U_g_size,
171-
user, &mat) );
172-
PetscCall( MatShellSetOperation(mat, MATOP_MULT,
173-
(void(*)(void))MatMult_Ceed) );
174-
PetscCall( MatShellSetVecType(mat, vec_type) );
175-
176-
KSP ksp;
177-
PetscCall( KSPCreate(comm, &ksp) );
178-
PetscCall( KSPSetOperators(ksp, mat, mat) );
179-
PetscCall( KSPSetFromOptions(ksp) );
180-
PetscCall( KSPSetUp(ksp) );
181-
PetscCall( KSPSolve(ksp, rhs, U_g) );
182-
//VecView(U_g, PETSC_VIEWER_STDOUT_WORLD);
183-
// ---------------------------------------------------------------------------
184-
// Compute pointwise L2 error
185-
// ---------------------------------------------------------------------------
186-
CeedScalar l2_error_u, l2_error_p;
187-
PetscCall( ComputeError(user, U_g, target,
188-
&l2_error_u, &l2_error_p) );
189147

190148
// ---------------------------------------------------------------------------
191-
// Output results
149+
// Solve PDE
192150
// ---------------------------------------------------------------------------
193-
KSPType ksp_type;
194-
KSPConvergedReason reason;
195-
PetscReal rnorm;
196-
PetscInt its;
197-
PetscCall( KSPGetType(ksp, &ksp_type) );
198-
PetscCall( KSPGetConvergedReason(ksp, &reason) );
199-
PetscCall( KSPGetIterationNumber(ksp, &its) );
200-
PetscCall( KSPGetResidualNorm(ksp, &rnorm) );
201-
PetscCall( PetscPrintf(comm,
202-
" KSP:\n"
203-
" KSP Type : %s\n"
204-
" KSP Convergence : %s\n"
205-
" Total KSP Iterations : %" PetscInt_FMT "\n"
206-
" Final rnorm : %e\n"
207-
" L2 Error of u and p : %e, %e\n",
208-
ksp_type, KSPConvergedReasons[reason], its,
209-
(double)rnorm, (double)l2_error_u,
210-
(double)l2_error_p) );
151+
Vec U_g;
152+
PDESolver(comm, dm, ceed, ceed_data, user, vec_type, target, rhs, &U_g);
153+
//VecView(U_g, PETSC_VIEWER_STDOUT_WORLD);
211154

212155
// ---------------------------------------------------------------------------
213156
// Save solution (paraview)
@@ -225,12 +168,8 @@ int main(int argc, char **argv) {
225168
// Free PETSc objects
226169
PetscCall( DMDestroy(&dm) );
227170
PetscCall( VecDestroy(&U_g) );
228-
PetscCall( VecDestroy(&U_loc) );
229171
PetscCall( VecDestroy(&rhs) );
230172
PetscCall( VecDestroy(&rhs_loc) );
231-
PetscCall( VecDestroy(&user->Y_loc) );
232-
PetscCall( MatDestroy(&mat) );
233-
PetscCall( KSPDestroy(&ksp) );
234173

235174
// -- Function list
236175
PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );

‎examples/Hdiv-mixed/main.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,5 +7,6 @@
77
#include "include/cl-options.h"
88
#include "include/problems.h"
99
#include "include/matops.h"
10+
#include "include/setup-solver.h"
1011

1112
#endif // MAIN_H

‎examples/Hdiv-mixed/src/matops.c‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -69,7 +69,7 @@ PetscErrorCode SNESFormResidual_Ceed(SNES snes, Vec X, Vec Y, void *ctx) {
6969
1.0, NULL, NULL, NULL) );
7070

7171
// libCEED for local action of residual evaluator
72-
PetscCall( ApplyLocalCeedOp(X, Y, user) );
72+
PetscCall( ApplyLocal_Ceed(user, X, Y) );
7373

7474
PetscFunctionReturn(0);
7575
};

‎examples/Hdiv-mixed/src/setup-libceed.c‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -146,8 +146,8 @@ PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm,
146146
// Set up libCEED on the fine grid for a given degree
147147
// -----------------------------------------------------------------------------
148148
PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
149-
ProblemData *problem_data, PetscInt U_g_size,
150-
PetscInt U_loc_size, CeedData ceed_data,
149+
ProblemData *problem_data,
150+
PetscInt rhs_loc_size, CeedData ceed_data,
151151
CeedVector rhs_ceed, CeedVector *target) {
152152
CeedInt P = app_ctx->degree + 1;
153153
// Number of quadratures in 1D, q_extra is set in cl-options.c
@@ -263,8 +263,8 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
263263
// Persistent libCEED vectors
264264
// ---------------------------------------------------------------------------
265265
// -- Operator action variables
266-
CeedVectorCreate(ceed, U_loc_size, &ceed_data->x_ceed);
267-
CeedVectorCreate(ceed, U_loc_size, &ceed_data->y_ceed);
266+
CeedVectorCreate(ceed, rhs_loc_size, &ceed_data->x_ceed);
267+
CeedVectorCreate(ceed, rhs_loc_size, &ceed_data->y_ceed);
268268

269269
// Local residual evaluator
270270
// ---------------------------------------------------------------------------

‎examples/Hdiv-mixed/src/setup-solver.c‎

Lines changed: 86 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,86 @@
1+
#include "../include/setup-solver.h"
2+
3+
// ---------------------------------------------------------------------------
4+
// Setup Solver
5+
// ---------------------------------------------------------------------------
6+
PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed,
7+
CeedData ceed_data, User user, VecType vec_type,
8+
CeedVector target, Vec rhs, Vec *U_g) {
9+
Vec U_loc;
10+
PetscInt U_l_size, U_g_size;
11+
12+
PetscFunctionBeginUser;
13+
14+
// Create global unknown solution U_g
15+
PetscCall( DMCreateGlobalVector(dm, U_g) );
16+
PetscCall( VecGetSize(*U_g, &U_g_size) );
17+
// Local size for matShell
18+
PetscCall( VecGetLocalSize(*U_g, &U_l_size) );
19+
// Create local unknown vector U_loc
20+
PetscCall( DMCreateLocalVector(dm, &U_loc) );
21+
22+
// ---------------------------------------------------------------------------
23+
// Setup Mat, KSP
24+
// ---------------------------------------------------------------------------
25+
user->comm = comm;
26+
user->dm = dm;
27+
user->X_loc = U_loc;
28+
PetscCall( VecDuplicate(U_loc, &user->Y_loc) );
29+
user->x_ceed = ceed_data->x_ceed;
30+
user->y_ceed = ceed_data->y_ceed;
31+
user->op_apply = ceed_data->op_residual;
32+
user->op_error = ceed_data->op_error;
33+
user->elem_restr_u = ceed_data->elem_restr_u;
34+
user->ceed = ceed;
35+
// Operator
36+
Mat mat;
37+
PetscCall( MatCreateShell(comm, U_l_size, U_l_size, U_g_size, U_g_size,
38+
user, &mat) );
39+
PetscCall( MatShellSetOperation(mat, MATOP_MULT,
40+
(void(*)(void))MatMult_Ceed) );
41+
PetscCall( MatShellSetVecType(mat, vec_type) );
42+
43+
KSP ksp;
44+
PetscCall( KSPCreate(comm, &ksp) );
45+
PetscCall( KSPSetOperators(ksp, mat, mat) );
46+
PetscCall( KSPSetFromOptions(ksp) );
47+
PetscCall( KSPSetUp(ksp) );
48+
PetscCall( KSPSolve(ksp, rhs, *U_g) );
49+
//VecView(*U_g, PETSC_VIEWER_STDOUT_WORLD);
50+
// ---------------------------------------------------------------------------
51+
// Compute pointwise L2 error
52+
// ---------------------------------------------------------------------------
53+
CeedScalar l2_error_u, l2_error_p;
54+
PetscCall( ComputeError(user, *U_g, target,
55+
&l2_error_u, &l2_error_p) );
56+
57+
// ---------------------------------------------------------------------------
58+
// Print output results
59+
// ---------------------------------------------------------------------------
60+
KSPType ksp_type;
61+
KSPConvergedReason reason;
62+
PetscReal rnorm;
63+
PetscInt its;
64+
PetscCall( KSPGetType(ksp, &ksp_type) );
65+
PetscCall( KSPGetConvergedReason(ksp, &reason) );
66+
PetscCall( KSPGetIterationNumber(ksp, &its) );
67+
PetscCall( KSPGetResidualNorm(ksp, &rnorm) );
68+
PetscCall( PetscPrintf(comm,
69+
" KSP:\n"
70+
" KSP Type : %s\n"
71+
" KSP Convergence : %s\n"
72+
" Total KSP Iterations : %" PetscInt_FMT "\n"
73+
" Final rnorm : %e\n"
74+
" L2 Error of u and p : %e, %e\n",
75+
ksp_type, KSPConvergedReasons[reason], its,
76+
(double)rnorm, (double)l2_error_u,
77+
(double)l2_error_p) );
78+
79+
// Free PETSc objects
80+
PetscCall( VecDestroy(&U_loc) );
81+
PetscCall( VecDestroy(&user->Y_loc) );
82+
PetscCall( MatDestroy(&mat) );
83+
PetscCall( KSPDestroy(&ksp) );
84+
85+
PetscFunctionReturn(0);
86+
};

0 commit comments

Comments
 (0)