Skip to content

Commit de6a9a6

Browse files
rezgarshakerirezgarshakeri
committed
Added visualization for Richards problem
1 parent f30f670 commit de6a9a6

File tree

10 files changed

+377
-287
lines changed

10 files changed

+377
-287
lines changed

‎examples/Hdiv-mixed/include/post-processing.h‎

Lines changed: 5 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -12,11 +12,10 @@ PetscErrorCode PrintOutput(Ceed ceed, AppCtx app_ctx, PetscBool has_ts,
1212
Vec U, CeedScalar l2_error_u,
1313
CeedScalar l2_error_p);
1414
PetscErrorCode SetupProjectVelocityCtx_Hdiv(MPI_Comm comm, DM dm, Ceed ceed,
15-
CeedData ceed_data, OperatorApplyContext ctx_post_Hdiv);
15+
CeedData ceed_data, OperatorApplyContext ctx_Hdiv);
1616
PetscErrorCode SetupProjectVelocityCtx_H1(MPI_Comm comm, DM dm_H1, Ceed ceed,
17-
CeedData ceed_data, OperatorApplyContext ctx_post_H1);
18-
PetscErrorCode ProjectVelocity(CeedData ceed_data,
19-
Vec U, VecType vec_type, Vec *U_H1,
20-
OperatorApplyContext ctx_post_Hdiv,
21-
OperatorApplyContext ctx_post_H1);
17+
CeedData ceed_data, VecType vec_type, OperatorApplyContext ctx_H1);
18+
PetscErrorCode ProjectVelocity(AppCtx app_ctx,
19+
Vec U, Vec *U_H1);
20+
PetscErrorCode CtxVecDestroy(AppCtx app_ctx);
2221
#endif // post_processing_h

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

Lines changed: 5 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -4,9 +4,11 @@
44
#include <ceed.h>
55
#include <petsc.h>
66

7+
#include "petscvec.h"
78
#include "structs.h"
89

910
PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
11+
VecType vec_type,
1012
OperatorApplyContext ctx_jacobian);
1113
PetscErrorCode SetupResidualOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
1214
OperatorApplyContext ctx_residual);
@@ -15,14 +17,9 @@ PetscErrorCode SetupErrorOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
1517
PetscErrorCode ApplyMatOp(Mat A, Vec X, Vec Y);
1618
PetscErrorCode SNESFormResidual(SNES snes, Vec X, Vec Y, void *ctx);
1719
PetscErrorCode SNESFormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx);
18-
PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data,
19-
VecType vec_type, SNES snes, KSP ksp, Vec *U);
20-
PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U,
20+
PetscErrorCode PDESolver(CeedData ceed_data, AppCtx app_ctx,
21+
SNES snes, KSP ksp, Vec *U);
22+
PetscErrorCode ComputeL2Error(CeedData ceed_data, AppCtx app_ctx, Vec U,
2123
CeedScalar *l2_error_u, CeedScalar *l2_error_p);
22-
PetscErrorCode PrintOutput(Ceed ceed, AppCtx app_ctx, PetscBool has_ts,
23-
CeedMemType mem_type_backend,
24-
TS ts, SNES snes, KSP ksp,
25-
Vec U, CeedScalar l2_error_u,
26-
CeedScalar l2_error_p);
2724

2825
#endif // setup_solvers_h

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

Lines changed: 4 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -6,11 +6,8 @@
66

77
#include "structs.h"
88

9-
PetscErrorCode CreateInitialConditions(CeedData ceed_data,
10-
Vec U, VecType vec_type,
11-
OperatorApplyContext ctx_initial_u0,
12-
OperatorApplyContext ctx_initial_p0,
13-
OperatorApplyContext ctx_residual_ut);
9+
PetscErrorCode CreateInitialConditions(CeedData ceed_data, AppCtx app_ctx,
10+
VecType vec_type, Vec U);
1411
PetscErrorCode SetupResidualOperatorCtx_Ut(MPI_Comm comm, DM dm, Ceed ceed,
1512
CeedData ceed_data, OperatorApplyContext ctx_residual_ut);
1613
PetscErrorCode SetupResidualOperatorCtx_U0(MPI_Comm comm, DM dm, Ceed ceed,
@@ -19,7 +16,7 @@ PetscErrorCode SetupResidualOperatorCtx_P0(MPI_Comm comm, DM dm, Ceed ceed,
1916
CeedData ceed_data, OperatorApplyContext ctx_initial_p0);
2017
PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
2118
void *ctx_residual_ut);
22-
PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx,
23-
Vec *U, TS *ts);
19+
PetscErrorCode TSSolveRichard(CeedData ceed_data, AppCtx app_ctx,
20+
TS ts, Vec *U);
2421

2522
#endif // setup_ts_h

‎examples/Hdiv-mixed/include/structs.h‎

Lines changed: 31 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -4,36 +4,20 @@
44
#include <ceed.h>
55
#include <petsc.h>
66

7-
// Application context from user command line options
8-
typedef struct AppCtx_ *AppCtx;
9-
struct AppCtx_ {
10-
MPI_Comm comm;
11-
// Degree of polynomial (1 only), extra quadrature pts
12-
PetscInt degree;
13-
PetscInt q_extra;
14-
// For applying traction BCs
15-
PetscInt bc_pressure_count;
16-
PetscInt bc_faces[16]; // face ID
17-
PetscScalar bc_pressure_value[16];
18-
// Problem type arguments
19-
PetscFunctionList problems;
20-
char problem_name[PETSC_MAX_PATH_LEN];
21-
CeedScalar t_final, t;
22-
PetscBool view_solution, quartic;
23-
};
24-
257
// PETSc operator contexts
268
typedef struct OperatorApplyContext_ *OperatorApplyContext;
279
struct OperatorApplyContext_ {
2810
MPI_Comm comm;
2911
Vec X_loc, Y_loc, X_t_loc;
30-
CeedVector x_ceed, y_ceed, x_t_ceed;
31-
CeedOperator op_apply;
12+
CeedVector x_ceed, y_ceed, x_t_ceed, x_coord, rhs_ceed_H1;
13+
CeedOperator op_apply, op_rhs_H1;
3214
DM dm;
3315
Ceed ceed;
3416
CeedScalar t, dt;
3517
CeedContextFieldLabel solution_time_label, final_time_label,
36-
timestep_label; ;
18+
timestep_label;
19+
CeedElemRestriction elem_restr_u_H1;
20+
VecType vec_type;
3721
};
3822

3923
// libCEED data struct
@@ -42,19 +26,39 @@ struct CeedData_ {
4226
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
4327
CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_U_i,
4428
elem_restr_p, elem_restr_p_i, elem_restr_u0,
45-
elem_restr_p0, elem_restr_u_post;
29+
elem_restr_p0, elem_restr_u_H1;
4630
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics_u, qf_ics_p,
47-
qf_rhs_u0, qf_rhs_p0, qf_post_rhs, qf_post_mass;
31+
qf_rhs_u0, qf_rhs_p0, qf_rhs_H1, qf_post_mass;
4832
CeedOperator op_residual, op_jacobian, op_error, op_ics_u, op_ics_p,
49-
op_rhs_u0, op_rhs_p0, op_post_rhs, op_post_mass;
33+
op_rhs_u0, op_rhs_p0, op_rhs_H1, op_post_mass;
5034
CeedVector x_ceed, y_ceed, x_coord, x_t_ceed, rhs_u0_ceed,
5135
u0_ceed, v0_ceed, rhs_p0_ceed, p0_ceed, q0_ceed,
52-
post_rhs_ceed, u_ceed, up_ceed, vp_ceed;
53-
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut,
54-
ctx_initial_u0, ctx_initial_p0, ctx_post_Hdiv, ctx_post_H1;
36+
rhs_ceed_H1, u_ceed, up_ceed, vp_ceed;
5537
CeedInt num_elem;
5638
};
5739

40+
// Application context from user command line options
41+
typedef struct AppCtx_ *AppCtx;
42+
struct AppCtx_ {
43+
MPI_Comm comm;
44+
// Degree of polynomial (1 only), extra quadrature pts
45+
PetscInt degree;
46+
PetscInt q_extra;
47+
// For applying traction BCs
48+
PetscInt bc_pressure_count;
49+
PetscInt bc_faces[16]; // face ID
50+
PetscScalar bc_pressure_value[16];
51+
// Problem type arguments
52+
PetscFunctionList problems;
53+
char problem_name[PETSC_MAX_PATH_LEN];
54+
CeedScalar t_final, t;
55+
PetscBool view_solution, quartic;
56+
char output_dir[PETSC_MAX_PATH_LEN];
57+
PetscInt output_freq;
58+
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut,
59+
ctx_initial_u0, ctx_initial_p0, ctx_Hdiv, ctx_H1;
60+
};
61+
5862
// Problem specific data
5963
typedef struct ProblemData_ *ProblemData;
6064
struct ProblemData_ {

‎examples/Hdiv-mixed/main.c‎

Lines changed: 77 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -73,18 +73,32 @@ int main(int argc, char **argv) {
7373
CeedData ceed_data;
7474
PetscCall( PetscCalloc1(1, &ceed_data) );
7575

76-
OperatorApplyContext ctx_residual_ut, ctx_initial_u0, ctx_initial_p0,
77-
ctx_post_Hdiv, ctx_post_H1;
76+
OperatorApplyContext ctx_jacobian, ctx_residual, ctx_residual_ut,
77+
ctx_initial_u0, ctx_initial_p0,
78+
ctx_error, ctx_Hdiv, ctx_H1;
79+
PetscCall( PetscCalloc1(1, &ctx_jacobian) );
80+
PetscCall( PetscCalloc1(1, &ctx_residual) );
7881
PetscCall( PetscCalloc1(1, &ctx_residual_ut) );
7982
PetscCall( PetscCalloc1(1, &ctx_initial_u0) );
8083
PetscCall( PetscCalloc1(1, &ctx_initial_p0) );
81-
PetscCall( PetscCalloc1(1, &ctx_post_Hdiv) );
82-
PetscCall( PetscCalloc1(1, &ctx_post_H1) );
83-
ceed_data->ctx_residual_ut = ctx_residual_ut;
84-
ceed_data->ctx_initial_u0 = ctx_initial_u0;
85-
ceed_data->ctx_initial_p0 = ctx_initial_p0;
86-
ceed_data->ctx_post_Hdiv = ctx_post_Hdiv;
87-
ceed_data->ctx_post_H1 = ctx_post_H1;
84+
PetscCall( PetscCalloc1(1, &ctx_error) );
85+
PetscCall( PetscCalloc1(1, &ctx_Hdiv) );
86+
PetscCall( PetscCalloc1(1, &ctx_H1) );
87+
// Context for Richards problem
88+
app_ctx->ctx_residual = ctx_residual;
89+
app_ctx->ctx_jacobian = ctx_jacobian;
90+
// Context for Richards problem
91+
app_ctx->ctx_residual_ut = ctx_residual_ut;
92+
// Context for initial velocity
93+
app_ctx->ctx_initial_u0 = ctx_initial_u0;
94+
// Context for initial pressure
95+
app_ctx->ctx_initial_p0 = ctx_initial_p0;
96+
// Context for MMS
97+
app_ctx->ctx_error = ctx_error;
98+
// Context for post-processing
99+
app_ctx->ctx_Hdiv = ctx_Hdiv;
100+
app_ctx->ctx_H1 = ctx_H1;
101+
88102
// ---------------------------------------------------------------------------
89103
// Process command line options
90104
// ---------------------------------------------------------------------------
@@ -112,23 +126,28 @@ int main(int argc, char **argv) {
112126
// Create DM
113127
// ---------------------------------------------------------------------------
114128
DM dm, dm_u0, dm_p0, dm_H1;
129+
// DM for mixed problem
115130
PetscCall( CreateDM(comm, vec_type, &dm) );
131+
// DM for projecting initial velocity to Hdiv space
116132
PetscCall( CreateDM(comm, vec_type, &dm_u0) );
133+
// DM for projecting initial pressure in L2
117134
PetscCall( CreateDM(comm, vec_type, &dm_p0) );
135+
// DM for projecting solution U into H1 space for PetscViewer
118136
PetscCall( CreateDM(comm, vec_type, &dm_H1) );
119137
// TODO: add mesh option
120138
// perturb to have smooth random mesh
121139
//PetscCall( PerturbVerticesSmooth(dm) );
122140

123141
// ---------------------------------------------------------------------------
124-
// Setup FE
142+
// Setup FE for H(div) mixed-problem and H1 projection in post-processing.c
125143
// ---------------------------------------------------------------------------
126-
SetupFEHdiv(comm, dm, dm_u0, dm_p0);
127-
SetupFEH1(problem_data, app_ctx, dm_H1);
144+
PetscCall( SetupFEHdiv(comm, dm, dm_u0, dm_p0) );
145+
PetscCall( SetupFEH1(problem_data, app_ctx, dm_H1) );
146+
128147
// ---------------------------------------------------------------------------
129-
// Create local Force vector
148+
// Create global unkown solution
130149
// ---------------------------------------------------------------------------
131-
Vec U; // U=[p,u], U0=u0
150+
Vec U; // U=[p,u]
132151
PetscCall( DMCreateGlobalVector(dm, &U) );
133152

134153
// ---------------------------------------------------------------------------
@@ -141,86 +160,89 @@ int main(int argc, char **argv) {
141160
//PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
142161
// bc_pressure) );
143162

163+
// ---------------------------------------------------------------------------
164+
// Setup context for projection problem; post-processing.c
165+
// ---------------------------------------------------------------------------
166+
PetscCall( SetupProjectVelocityCtx_Hdiv(comm, dm, ceed, ceed_data,
167+
app_ctx->ctx_Hdiv) );
168+
PetscCall( SetupProjectVelocityCtx_H1(comm, dm_H1, ceed, ceed_data,
169+
vec_type, app_ctx->ctx_H1) );
170+
144171
// ---------------------------------------------------------------------------
145172
// Setup TSSolve for Richard problem
146173
// ---------------------------------------------------------------------------
147174
TS ts;
148175
if (problem_data->has_ts) {
149176
// ---------------------------------------------------------------------------
150-
// Create global initial conditions
177+
// Setup context for initial conditions
151178
// ---------------------------------------------------------------------------
152-
SetupResidualOperatorCtx_U0(comm, dm_u0, ceed, ceed_data,
153-
ceed_data->ctx_initial_u0);
154-
SetupResidualOperatorCtx_P0(comm, dm_p0, ceed, ceed_data,
155-
ceed_data->ctx_initial_p0);
156-
SetupResidualOperatorCtx_Ut(comm, dm, ceed, ceed_data,
157-
ceed_data->ctx_residual_ut);
158-
CreateInitialConditions(ceed_data, U, vec_type,
159-
ceed_data->ctx_initial_u0,
160-
ceed_data->ctx_initial_p0,
161-
ceed_data->ctx_residual_ut);
179+
PetscCall( SetupResidualOperatorCtx_U0(comm, dm_u0, ceed, ceed_data,
180+
app_ctx->ctx_initial_u0) );
181+
PetscCall( SetupResidualOperatorCtx_P0(comm, dm_p0, ceed, ceed_data,
182+
app_ctx->ctx_initial_p0) );
183+
PetscCall( SetupResidualOperatorCtx_Ut(comm, dm, ceed, ceed_data,
184+
app_ctx->ctx_residual_ut) );
185+
PetscCall( CreateInitialConditions(ceed_data, app_ctx, vec_type, U) );
162186
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
163187
// Solve Richards problem
164-
PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_loc) );
165-
PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_t_loc) );
166-
PetscCall( TSSolveRichard(dm, ceed_data, app_ctx,
167-
&U, &ts) );
188+
PetscCall( TSCreate(comm, &ts) );
189+
PetscCall( VecZeroEntries(app_ctx->ctx_residual_ut->X_loc) );
190+
PetscCall( VecZeroEntries(app_ctx->ctx_residual_ut->X_t_loc) );
191+
PetscCall( TSSolveRichard(ceed_data, app_ctx, ts, &U) );
168192
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
169193
}
170194

195+
// ---------------------------------------------------------------------------
196+
// Setup SNES for Darcy problem
197+
// ---------------------------------------------------------------------------
171198
SNES snes;
172199
KSP ksp;
173200
if (!problem_data->has_ts) {
174-
// ---------------------------------------------------------------------------
175-
// Setup SNES for Darcy problem
176-
// ---------------------------------------------------------------------------
201+
PetscCall( SetupJacobianOperatorCtx(dm, ceed, ceed_data, vec_type,
202+
app_ctx->ctx_jacobian) );
203+
PetscCall( SetupResidualOperatorCtx(dm, ceed, ceed_data,
204+
app_ctx->ctx_residual) );
177205
// Create SNES
178206
PetscCall( SNESCreate(comm, &snes) );
179207
PetscCall( SNESGetKSP(snes, &ksp) );
180-
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
208+
PetscCall( PDESolver(ceed_data, app_ctx, snes, ksp, &U) );
181209
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
182210
}
183211

184212
// ---------------------------------------------------------------------------
185213
// Compute L2 error of mms problem
186214
// ---------------------------------------------------------------------------
215+
PetscCall( SetupErrorOperatorCtx(dm, ceed, ceed_data, app_ctx->ctx_error) );
187216
CeedScalar l2_error_u, l2_error_p;
188-
PetscCall( ComputeL2Error(dm, ceed, ceed_data, U, &l2_error_u,
189-
&l2_error_p) );
217+
PetscCall( ComputeL2Error(ceed_data, app_ctx, U,
218+
&l2_error_u, &l2_error_p) );
190219

191220
// ---------------------------------------------------------------------------
192-
// Print output results
221+
// Print solver iterations and final norms
193222
// ---------------------------------------------------------------------------
194223
PetscCall( PrintOutput(ceed, app_ctx, problem_data->has_ts, mem_type_backend,
195224
ts, snes, ksp, U, l2_error_u, l2_error_p) );
196225

197226
// ---------------------------------------------------------------------------
198227
// Save solution (paraview)
199228
// ---------------------------------------------------------------------------
200-
Vec U_H1;
201-
PetscCall( DMCreateGlobalVector(dm_H1, &U_H1) );
202-
PetscCall( VecZeroEntries(U_H1) );
203229
if (app_ctx->view_solution) {
204230
PetscViewer viewer_p;
205-
PetscCall( PetscViewerVTKOpen(comm,"solution_p.vtu",FILE_MODE_WRITE,
231+
PetscCall( PetscViewerVTKOpen(comm,"darcy_pressure.vtu",FILE_MODE_WRITE,
206232
&viewer_p) );
207233
PetscCall( VecView(U, viewer_p) );
208234
PetscCall( PetscViewerDestroy(&viewer_p) );
209235

210-
SetupProjectVelocityCtx_Hdiv(comm, dm, ceed, ceed_data,
211-
ceed_data->ctx_post_Hdiv);
212-
SetupProjectVelocityCtx_H1(comm, dm_H1, ceed, ceed_data,
213-
ceed_data->ctx_post_H1);
214-
215-
ProjectVelocity(ceed_data, U, vec_type, &U_H1,
216-
ceed_data->ctx_post_Hdiv,
217-
ceed_data->ctx_post_H1);
236+
Vec U_H1; // velocity in H1 space for post-processing
237+
PetscCall( DMCreateGlobalVector(dm_H1, &U_H1) );
238+
PetscCall( ProjectVelocity(app_ctx, U, &U_H1) );
218239

219240
PetscViewer viewer_u;
220-
PetscCall( PetscViewerVTKOpen(comm,"solution_u.vtu",FILE_MODE_WRITE,
241+
PetscCall( PetscViewerVTKOpen(comm,"darcy_velocity.vtu",FILE_MODE_WRITE,
221242
&viewer_u) );
222243
PetscCall( VecView(U_H1, viewer_u) );
223244
PetscCall( PetscViewerDestroy(&viewer_u) );
245+
PetscCall( VecDestroy(&U_H1) );
224246
}
225247
// ---------------------------------------------------------------------------
226248
// Free objects
@@ -232,10 +254,7 @@ int main(int argc, char **argv) {
232254
PetscCall( DMDestroy(&dm_p0) );
233255
PetscCall( DMDestroy(&dm_H1) );
234256
PetscCall( VecDestroy(&U) );
235-
PetscCall( VecDestroy(&U_H1) );
236-
PetscCall( VecDestroy(&ceed_data->ctx_residual_ut->X_loc) );
237-
PetscCall( VecDestroy(&ceed_data->ctx_residual_ut->X_t_loc) );
238-
PetscCall( VecDestroy(&ceed_data->ctx_residual_ut->Y_loc) );
257+
PetscCall( CtxVecDestroy(app_ctx) );
239258
if (problem_data->has_ts) {
240259
PetscCall( TSDestroy(&ts) );
241260
} else {
@@ -252,8 +271,11 @@ int main(int argc, char **argv) {
252271
PetscCall( PetscFree(ctx_initial_u0) );
253272
PetscCall( PetscFree(ctx_initial_p0) );
254273
PetscCall( PetscFree(ctx_residual_ut) );
255-
PetscCall( PetscFree(ctx_post_H1) );
256-
PetscCall( PetscFree(ctx_post_Hdiv) );
274+
PetscCall( PetscFree(ctx_residual) );
275+
PetscCall( PetscFree(ctx_jacobian) );
276+
PetscCall( PetscFree(ctx_error) );
277+
PetscCall( PetscFree(ctx_H1) );
278+
PetscCall( PetscFree(ctx_Hdiv) );
257279

258280
// Free libCEED objects
259281
//CeedVectorDestroy(&bc_pressure);

0 commit comments

Comments
 (0)