Skip to content

Commit 5e6a9fe

Browse files
rezgarshakerirezgarshakeri
committed
Setup Ifunctiona and TSSolve for Richards problem
1 parent 833e410 commit 5e6a9fe

File tree

13 files changed

+445
-388
lines changed

13 files changed

+445
-388
lines changed

‎examples/Hdiv-mixed/conv_test_result.csv‎

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
run,mesh_res,error_u,error_p
22
0,2,9.13818,0.09469
3-
1,3,4.60103,0.05186
4-
2,4,2.75883,0.03199
5-
3,5,1.81813,0.02132
6-
4,6,1.27780,0.01505
7-
5,7,0.93759,0.01108
8-
6,8,0.71073,0.00842
9-
7,9,0.55232,0.00655
10-
8,10,0.43767,0.00520
11-
9,11,0.35229,0.00418
12-
10,12,0.28723,0.00340
3+
1,3,4.44275,0.05113
4+
2,4,2.55202,0.03052
5+
3,5,1.63121,0.01989
6+
4,6,1.11886,0.01381
7+
5,7,0.80594,0.01003
8+
6,8,0.60143,0.00754
9+
7,9,0.46079,0.00581
10+
8,10,0.36019,0.00456
11+
9,11,0.28599,0.00364
12+
10,12,0.22996,0.00293

‎examples/Hdiv-mixed/convrate_mixed.png‎

573 Bytes
Loading

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

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@
66
// Convert PETSc MemType to libCEED MemType
77
CeedMemType MemTypeP2C(PetscMemType mtype);
88
// Destroy libCEED objects
9-
PetscErrorCode CeedDataDestroy(CeedData ceed_data);
9+
PetscErrorCode CeedDataDestroy(CeedData ceed_data, ProblemData problem_data);
1010
// Utility function - essential BC dofs are encoded in closure indices as -(i+1)
1111
PetscInt Involute(PetscInt i);
1212
// Utility function to create local CEED restriction from DMPlex
@@ -18,6 +18,6 @@ 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,
22-
PetscInt U_loc_size, CeedData ceed_data);
21+
OperatorApplyContext ctx_residual_ut,
22+
ProblemData problem_data, CeedData ceed_data);
2323
#endif // setuplibceed_h

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

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -20,9 +20,10 @@ PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data,
2020
PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U,
2121
CeedScalar *l2_error_u, CeedScalar *l2_error_p);
2222
PetscErrorCode PrintOutput(Ceed ceed,
23-
CeedMemType mem_type_backend,
23+
CeedMemType mem_type_backend, TS ts,
2424
SNES snes, KSP ksp,
2525
Vec U, CeedScalar l2_error_u,
26-
CeedScalar l2_error_p, AppCtx app_ctx);
26+
CeedScalar l2_error_p, AppCtx app_ctx,
27+
PetscBool has_ts);
2728

2829
#endif // setup_solvers_h

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

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -6,12 +6,14 @@
66

77
#include "structs.h"
88

9-
PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *U0);
9+
PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *Y0,
10+
OperatorApplyContext ctx_residual_ut);
1011
PetscErrorCode SetupResidualOperatorCtx_Ut(DM dm, Ceed ceed, CeedData ceed_data,
1112
OperatorApplyContext ctx_residual_ut);
1213
PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
1314
void *ctx_residual_ut);
14-
PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx,
15-
Vec *U, PetscScalar *f_time, TS *ts);
15+
PetscErrorCode TSSolveRichard(DM dm, Ceed ceed, CeedData ceed_data,
16+
AppCtx app_ctx, OperatorApplyContext ctx_residual_ut,
17+
Vec *U, TS *ts);
1618

1719
#endif // setup_ts_h

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

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
#ifndef structs_h
22
#define structs_h
33

4+
#include "ceed/ceed-f64.h"
5+
#include "petscsystypes.h"
46
#include <ceed.h>
57
#include <petsc.h>
68

@@ -18,7 +20,7 @@ struct AppCtx_ {
1820
// Problem type arguments
1921
PetscFunctionList problems;
2022
char problem_name[PETSC_MAX_PATH_LEN];
21-
CeedContextFieldLabel solution_time_label;
23+
CeedScalar t_final;
2224
};
2325

2426
// 2) richard
@@ -70,9 +72,11 @@ struct OperatorApplyContext_ {
7072
MPI_Comm comm;
7173
Vec X_loc, Y_loc, X_t_loc;
7274
CeedVector x_ceed, y_ceed, x_t_ceed;
73-
CeedOperator op_apply;
75+
CeedOperator op_apply, op_true;
7476
DM dm;
7577
Ceed ceed;
78+
CeedScalar t;
79+
CeedContextFieldLabel solution_time_label, final_time_label;
7680
};
7781

7882
// libCEED data struct
@@ -81,10 +85,10 @@ struct CeedData_ {
8185
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
8286
CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_U_i,
8387
elem_restr_p, elem_restr_p_i;
84-
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics;
85-
CeedOperator op_residual, op_jacobian, op_error, op_ics;
86-
CeedVector x_ceed, y_ceed, x_coord, U0_ceed, x_t_ceed;
87-
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut;
88+
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics, qf_true;
89+
CeedOperator op_residual, op_jacobian, op_error, op_ics, op_true;
90+
CeedVector x_ceed, y_ceed, x_coord, x_t_ceed;
91+
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error;
8892
};
8993

9094
// Problem specific data

‎examples/Hdiv-mixed/main.c‎

Lines changed: 48 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -28,6 +28,10 @@
2828
// ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
2929
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
3030
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
31+
// ./main -pc_type svd -problem richard2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
32+
33+
#include "ceed/ceed.h"
34+
#include <stdio.h>
3135
const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
3236

3337
#include "main.h"
@@ -76,6 +80,9 @@ int main(int argc, char **argv) {
7680
Physics phys_ctx;
7781
PetscCall( PetscCalloc1(1, &phys_ctx) );
7882

83+
OperatorApplyContext ctx_residual_ut;
84+
PetscCall( PetscCalloc1(1, &ctx_residual_ut) );
85+
7986
// ---------------------------------------------------------------------------
8087
// Process command line options
8188
// ---------------------------------------------------------------------------
@@ -113,65 +120,66 @@ int main(int argc, char **argv) {
113120
// ---------------------------------------------------------------------------
114121
SetupFE(comm, dm);
115122

116-
// ---------------------------------------------------------------------------
117-
// Create local Force vector
118-
// ---------------------------------------------------------------------------
119-
Vec U_loc;
120-
PetscInt U_loc_size;
121-
//CeedVector bc_pressure;
122-
PetscCall( DMCreateLocalVector(dm, &U_loc) );
123-
// Local size for libCEED
124-
PetscCall( VecGetSize(U_loc, &U_loc_size) );
125-
126123
// ---------------------------------------------------------------------------
127124
// Setup libCEED
128125
// ---------------------------------------------------------------------------
129126
// -- Set up libCEED objects
130-
PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
131-
U_loc_size, ceed_data) );
127+
PetscCall( SetupLibceed(dm, ceed, app_ctx, ctx_residual_ut,
128+
problem_data, ceed_data) );
132129
//CeedVectorView(force_ceed, "%12.8f", stdout);
133130
//PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
134131
// bc_pressure) );
135132

136133

134+
// ---------------------------------------------------------------------------
135+
// Create Global Solution
136+
// ---------------------------------------------------------------------------
137+
Vec U; // U = [p,u]
138+
PetscCall( DMCreateGlobalVector(dm, &U) );
139+
137140
// ---------------------------------------------------------------------------
138141
// Setup TSSolve for Richard problem
139142
// ---------------------------------------------------------------------------
143+
TS ts;
144+
SNES snes;
145+
KSP ksp;
140146
if (problem_data->has_ts) {
141147
// ---------------------------------------------------------------------------
142148
// Create global initial conditions
143149
// ---------------------------------------------------------------------------
144-
Vec U0;
145-
CreateInitialConditions(dm, ceed_data, &U0);
146-
VecView(U0, PETSC_VIEWER_STDOUT_WORLD);
147-
PetscCall( VecDestroy(&U0) );
150+
SetupResidualOperatorCtx_Ut(dm, ceed, ceed_data, ctx_residual_ut);
151+
152+
//SetupResidualOperatorCtx_U0(dm, ceed, ceed_data,
153+
// ctx_initial);
154+
CreateInitialConditions(dm, ceed_data, &U, ctx_residual_ut);
155+
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
156+
PetscCall( VecZeroEntries(ctx_residual_ut->X_t_loc) );
157+
PetscCall( TSSolveRichard(dm, ceed, ceed_data, app_ctx, ctx_residual_ut,
158+
&U, &ts) );
148159
}
149160

150-
// ---------------------------------------------------------------------------
151-
// Solve PDE
152-
// ---------------------------------------------------------------------------
153-
// Create SNES
154-
SNES snes;
155-
KSP ksp;
156-
Vec U;
157-
PetscCall( SNESCreate(comm, &snes) );
158-
PetscCall( SNESGetKSP(snes, &ksp) );
159-
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
160-
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
161+
if (!problem_data->has_ts) {
162+
// ---------------------------------------------------------------------------
163+
// Solve PDE
164+
// ---------------------------------------------------------------------------
165+
// Create SNES
166+
PetscCall( SNESCreate(comm, &snes) );
167+
PetscCall( SNESGetKSP(snes, &ksp) );
168+
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
169+
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
170+
}
161171

162172
// ---------------------------------------------------------------------------
163173
// Compute L2 error of mms problem
164174
// ---------------------------------------------------------------------------
165175
CeedScalar l2_error_u, l2_error_p;
166176
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u,
167177
&l2_error_p) );
168-
169178
// ---------------------------------------------------------------------------
170179
// Print output results
171180
// ---------------------------------------------------------------------------
172-
PetscCall( PrintOutput(ceed, mem_type_backend,
173-
snes, ksp, U, l2_error_u, l2_error_p, app_ctx) );
174-
181+
PetscCall( PrintOutput(ceed, mem_type_backend, ts,
182+
snes, ksp, U, l2_error_u, l2_error_p, app_ctx, problem_data->has_ts) );
175183
// ---------------------------------------------------------------------------
176184
// Save solution (paraview)
177185
// ---------------------------------------------------------------------------
@@ -188,20 +196,23 @@ int main(int argc, char **argv) {
188196
// Free PETSc objects
189197
PetscCall( DMDestroy(&dm) );
190198
PetscCall( VecDestroy(&U) );
191-
PetscCall( VecDestroy(&U_loc) );
192-
PetscCall( SNESDestroy(&snes) );
199+
if (problem_data->has_ts) {
200+
PetscCall( TSDestroy(&ts) );
201+
} else {
202+
PetscCall( SNESDestroy(&snes) );
203+
}
193204

194205
// -- Function list
195206
PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );
196207

208+
// Free libCEED objects
209+
//CeedVectorDestroy(&bc_pressure);
210+
PetscCall( CeedDataDestroy(ceed_data, problem_data) );
197211
// -- Structs
198212
PetscCall( PetscFree(app_ctx) );
199213
PetscCall( PetscFree(problem_data) );
200214
PetscCall( PetscFree(phys_ctx) );
201-
202-
// Free libCEED objects
203-
//CeedVectorDestroy(&bc_pressure);
204-
PetscCall( CeedDataDestroy(ceed_data) );
215+
PetscCall( PetscFree(ctx_residual_ut) );
205216
CeedDestroy(&ceed);
206217

207218
return PetscFinalize();

‎examples/Hdiv-mixed/problems/richard2d.c‎

Lines changed: 16 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,7 @@
2020
#include "../include/register-problem.h"
2121
#include "../qfunctions/richard-system2d.h"
2222
#include "../qfunctions/richard-true2d.h"
23+
#include "../qfunctions/darcy-error2d.h"
2324
#include "../qfunctions/pressure-boundary2d.h"
2425
#include "petscsystypes.h"
2526

@@ -43,12 +44,12 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
4344
problem_data->ics_loc = RichardICs2D_loc;
4445
problem_data->true_solution = RichardTrue2D;
4546
problem_data->true_solution_loc = RichardTrue2D_loc;
46-
//problem_data->residual = RichardSystem2D;
47-
//problem_data->residual_loc = RichardSystem2D_loc;
48-
//problem_data->jacobian = JacobianRichardSystem2D;
49-
//problem_data->jacobian_loc = JacobianRichardSystem2D_loc;
50-
//problem_data->error = DarcyError2D;
51-
//problem_data->error_loc = DarcyError2D_loc;
47+
problem_data->residual = RichardSystem2D;
48+
problem_data->residual_loc = RichardSystem2D_loc;
49+
problem_data->jacobian = JacobianRichardSystem2D;
50+
problem_data->jacobian_loc = JacobianRichardSystem2D_loc;
51+
problem_data->error = DarcyError2D;
52+
problem_data->error_loc = DarcyError2D_loc;
5253
problem_data->bc_pressure = BCPressure2D;
5354
problem_data->bc_pressure_loc = BCPressure2D_loc;
5455
problem_data->has_ts = PETSC_TRUE;
@@ -72,6 +73,11 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
7273
rho_a0, &rho_a0, NULL));
7374
PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
7475
beta, &beta, NULL));
76+
PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
77+
beta, &beta, NULL));
78+
app_ctx->t_final = 1.;
79+
PetscCall( PetscOptionsScalar("-t_final", "End time", NULL,
80+
app_ctx->t_final, &app_ctx->t_final, NULL));
7581
PetscOptionsEnd();
7682

7783
richard_ctx->kappa = kappa;
@@ -82,11 +88,15 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
8288
richard_ctx->g = g;
8389
richard_ctx->p0 = p0;
8490
richard_ctx->gamma = 5.;
91+
richard_ctx->t = 0.;
92+
richard_ctx->t_final = app_ctx->t_final;
8593
CeedQFunctionContextCreate(ceed, &richard_context);
8694
CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
8795
sizeof(*richard_ctx), richard_ctx);
8896
CeedQFunctionContextRegisterDouble(richard_context, "time",
8997
offsetof(struct RICHARDContext_, t), 1, "current solver time");
98+
CeedQFunctionContextRegisterDouble(richard_context, "final_time",
99+
offsetof(struct RICHARDContext_, t_final), 1, "final time");
90100
problem_data->qfunction_context = richard_context;
91101
CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
92102
FreeContextPetsc);

0 commit comments

Comments
 (0)