Skip to content

Commit 3a8fc58

Browse files
rezgarshakerirezgarshakeri
committed
Create Initial conditions for Richards problem
1 parent 833e410 commit 3a8fc58

25 files changed

+1029
-405
lines changed

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,6 @@
1010
// ---------------------------------------------------------------------------
1111
// Setup FE
1212
// ---------------------------------------------------------------------------
13-
PetscErrorCode SetupFE(MPI_Comm comm, DM dm);
13+
PetscErrorCode SetupFE(MPI_Comm comm, DM dm, DM dm_u0, DM dm_p0);
1414

1515
#endif // setupfe_h

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

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -6,18 +6,21 @@
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
1313
PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm,
1414
CeedInt height, DMLabel domain_label, CeedInt value,
1515
CeedElemRestriction *elem_restr);
1616
// Utility function to create local CEED Oriented restriction from DMPlex
17-
PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P,
18-
CeedElemRestriction *elem_restr_oriented, CeedElemRestriction *elem_restr);
17+
PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, DM dm_u0,
18+
DM dm_p0, CeedInt P,
19+
CeedElemRestriction *elem_restr_u, CeedElemRestriction *elem_restr_p,
20+
CeedElemRestriction *elem_restr_u0, CeedElemRestriction *elem_restr_p0);
1921
// Set up libCEED for a given degree
20-
PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
22+
PetscErrorCode SetupLibceed(DM dm, DM dm_u0, DM dm_p0, Ceed ceed,
23+
AppCtx app_ctx,
2124
ProblemData problem_data,
22-
PetscInt U_loc_size, CeedData ceed_data);
25+
CeedData ceed_data);
2326
#endif // setuplibceed_h

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -12,7 +12,7 @@ PetscErrorCode SetupResidualOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
1212
OperatorApplyContext ctx_residual);
1313
PetscErrorCode SetupErrorOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
1414
OperatorApplyContext ctx_error);
15-
PetscErrorCode ApplyJacobian(Mat A, Vec X, Vec Y);
15+
PetscErrorCode ApplyMatOp(Mat A, Vec X, Vec Y);
1616
PetscErrorCode SNESFormResidual(SNES snes, Vec X, Vec Y, void *ctx);
1717
PetscErrorCode SNESFormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx);
1818
PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data,

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

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,9 +6,19 @@
66

77
#include "structs.h"
88

9-
PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *U0);
9+
PetscErrorCode CreateInitialConditions(DM dm, DM dm_u0, DM dm_p0,
10+
CeedData ceed_data,
11+
Vec U, VecType vec_type,
12+
OperatorApplyContext ctx_initial_u0,
13+
OperatorApplyContext ctx_initial_p0);
1014
PetscErrorCode SetupResidualOperatorCtx_Ut(DM dm, Ceed ceed, CeedData ceed_data,
1115
OperatorApplyContext ctx_residual_ut);
16+
PetscErrorCode SetupResidualOperatorCtx_U0(MPI_Comm comm, DM dm, Ceed ceed,
17+
CeedData ceed_data,
18+
OperatorApplyContext ctx_initial_u0);
19+
PetscErrorCode SetupResidualOperatorCtx_P0(MPI_Comm comm, DM dm, Ceed ceed,
20+
CeedData ceed_data,
21+
OperatorApplyContext ctx_initial_p0);
1222
PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
1323
void *ctx_residual_ut);
1424
PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx,

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

Lines changed: 18 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -19,49 +19,7 @@ struct AppCtx_ {
1919
PetscFunctionList problems;
2020
char problem_name[PETSC_MAX_PATH_LEN];
2121
CeedContextFieldLabel solution_time_label;
22-
};
23-
24-
// 2) richard
25-
// We have 3 experiment parameters as described in Table 1:P1, P2, P3
26-
// Matthew Farthing, Christopher Kees, Cass Miller (2003)
27-
// https://www.sciencedirect.com/science/article/pii/S0309170802001872
28-
29-
#ifndef PHYSICS_RICHARDP2_STRUCT
30-
#define PHYSICS_RICHARDP2_STRUCT
31-
typedef struct RICHARDP2Context_ *RICHARDP2Context;
32-
struct RICHARDP2Context_ {
33-
CeedScalar K_star;
34-
CeedScalar theta_s;
35-
CeedScalar theta_r;
36-
CeedScalar alpha_v;
37-
CeedScalar n_v;
38-
CeedScalar m_v;
39-
CeedScalar m_r;
40-
CeedScalar rho_0;
41-
CeedScalar beta;
42-
};
43-
#endif
44-
45-
#ifndef PHYSICS_RICHARDP3_STRUCT
46-
#define PHYSICS_RICHARDP3_STRUCT
47-
typedef struct RICHARDP3Context_ *RICHARDP3Context;
48-
struct RICHARDP3Context_ {
49-
CeedScalar K_star;
50-
CeedScalar theta_s;
51-
CeedScalar theta_r;
52-
CeedScalar alpha_star_v;
53-
CeedScalar n_v;
54-
CeedScalar m_v;
55-
CeedScalar rho_0;
56-
CeedScalar beta;
57-
};
58-
#endif
59-
60-
// Struct that contains all enums and structs used for the physics of all problems
61-
typedef struct Physics_ *Physics;
62-
struct Physics_ {
63-
RICHARDP2Context richard_p2_ctx;
64-
RICHARDP3Context richard_p3_ctx;
22+
CeedScalar t_final;
6523
};
6624

6725
// PETSc operator contexts
@@ -73,27 +31,36 @@ struct OperatorApplyContext_ {
7331
CeedOperator op_apply;
7432
DM dm;
7533
Ceed ceed;
34+
CeedScalar t;
35+
CeedContextFieldLabel solution_time_label, final_time_label;
7636
};
7737

7838
// libCEED data struct
7939
typedef struct CeedData_ *CeedData;
8040
struct CeedData_ {
8141
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
8242
CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_U_i,
83-
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;
43+
elem_restr_p, elem_restr_p_i, elem_restr_u0,
44+
elem_restr_p0;
45+
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics_u, qf_ics_p,
46+
qf_rhs_u0, qf_rhs_p0;
47+
CeedOperator op_residual, op_jacobian, op_error, op_ics_u, op_ics_p,
48+
op_rhs_u0, op_rhs_p0;
49+
CeedVector x_ceed, y_ceed, x_coord, x_t_ceed, rhs_u0_ceed,
50+
u0_ceed, v0_ceed, rhs_p0_ceed, p0_ceed, q0_ceed;
51+
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut,
52+
ctx_initial_u0, ctx_initial_p0;
53+
CeedInt num_elem;
8854
};
8955

9056
// Problem specific data
9157
typedef struct ProblemData_ *ProblemData;
9258
struct ProblemData_ {
93-
CeedQFunctionUser true_solution, residual, jacobian, error, ics,
94-
bc_pressure;
59+
CeedQFunctionUser true_solution, residual, jacobian, error, ics_u, ics_p,
60+
bc_pressure, rhs_u0, rhs_p0;
9561
const char *true_solution_loc, *residual_loc, *jacobian_loc,
96-
*error_loc, *bc_pressure_loc, *ics_loc;
62+
*error_loc, *bc_pressure_loc, *ics_u_loc, *ics_p_loc, *rhs_u0_loc,
63+
*rhs_p0_loc;
9764
CeedQuadMode quadrature_mode;
9865
CeedInt elem_node, dim, q_data_size_face;
9966
CeedQFunctionContext qfunction_context;

‎examples/Hdiv-mixed/main.c

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

76-
Physics phys_ctx;
77-
PetscCall( PetscCalloc1(1, &phys_ctx) );
78-
76+
OperatorApplyContext ctx_residual_ut, ctx_initial_u0, ctx_initial_p0;
77+
PetscCall( PetscCalloc1(1, &ctx_residual_ut) );
78+
PetscCall( PetscCalloc1(1, &ctx_initial_u0) );
79+
PetscCall( PetscCalloc1(1, &ctx_initial_p0) );
80+
ceed_data->ctx_residual_ut = ctx_residual_ut;
81+
ceed_data->ctx_initial_u0 = ctx_initial_u0;
82+
ceed_data->ctx_initial_p0 = ctx_initial_p0;
7983
// ---------------------------------------------------------------------------
8084
// Process command line options
8185
// ---------------------------------------------------------------------------
@@ -102,33 +106,32 @@ int main(int argc, char **argv) {
102106
// ---------------------------------------------------------------------------
103107
// Create DM
104108
// ---------------------------------------------------------------------------
105-
DM dm;
109+
DM dm, dm_u0, dm_p0;
106110
PetscCall( CreateDM(comm, vec_type, &dm) );
111+
PetscCall( CreateDM(comm, vec_type, &dm_u0) );
112+
PetscCall( CreateDM(comm, vec_type, &dm_p0) );
107113
// TODO: add mesh option
108114
// perturb to have smooth random mesh
109-
// PetscCall( PerturbVerticesSmooth(dm) );
115+
//PetscCall( PerturbVerticesSmooth(dm) );
110116

111117
// ---------------------------------------------------------------------------
112118
// Setup FE
113119
// ---------------------------------------------------------------------------
114-
SetupFE(comm, dm);
120+
SetupFE(comm, dm, dm_u0, dm_p0);
115121

116122
// ---------------------------------------------------------------------------
117123
// Create local Force vector
118124
// ---------------------------------------------------------------------------
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+
Vec U; // U=[p,u], U0=u0
126+
PetscCall( DMCreateGlobalVector(dm, &U) );
127+
PetscCall( VecZeroEntries(U) );
125128

126129
// ---------------------------------------------------------------------------
127130
// Setup libCEED
128131
// ---------------------------------------------------------------------------
129132
// -- Set up libCEED objects
130-
PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
131-
U_loc_size, ceed_data) );
133+
PetscCall( SetupLibceed(dm, dm_u0, dm_p0, ceed, app_ctx, problem_data,
134+
ceed_data) );
132135
//CeedVectorView(force_ceed, "%12.8f", stdout);
133136
//PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
134137
// bc_pressure) );
@@ -141,67 +144,76 @@ int main(int argc, char **argv) {
141144
// ---------------------------------------------------------------------------
142145
// Create global initial conditions
143146
// ---------------------------------------------------------------------------
144-
Vec U0;
145-
CreateInitialConditions(dm, ceed_data, &U0);
146-
VecView(U0, PETSC_VIEWER_STDOUT_WORLD);
147-
PetscCall( VecDestroy(&U0) );
147+
SetupResidualOperatorCtx_U0(comm, dm_u0, ceed, ceed_data,
148+
ceed_data->ctx_initial_u0);
149+
SetupResidualOperatorCtx_P0(comm, dm_u0, ceed, ceed_data,
150+
ceed_data->ctx_initial_p0);
151+
CreateInitialConditions(dm, dm_u0, dm_p0, ceed_data, U, vec_type,
152+
ceed_data->ctx_initial_u0,
153+
ceed_data->ctx_initial_p0);
154+
VecView(U, PETSC_VIEWER_STDOUT_WORLD);
148155
}
149156

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);
157+
if (!problem_data->has_ts) {
158+
// ---------------------------------------------------------------------------
159+
// Solve PDE
160+
// ---------------------------------------------------------------------------
161+
// Create SNES
162+
SNES snes;
163+
KSP ksp;
164+
PetscCall( SNESCreate(comm, &snes) );
165+
PetscCall( SNESGetKSP(snes, &ksp) );
166+
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
167+
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
161168

162-
// ---------------------------------------------------------------------------
163-
// Compute L2 error of mms problem
164-
// ---------------------------------------------------------------------------
165-
CeedScalar l2_error_u, l2_error_p;
166-
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u,
167-
&l2_error_p) );
169+
// ---------------------------------------------------------------------------
170+
// Compute L2 error of mms problem
171+
// ---------------------------------------------------------------------------
172+
CeedScalar l2_error_u, l2_error_p;
173+
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u,
174+
&l2_error_p) );
168175

169-
// ---------------------------------------------------------------------------
170-
// Print output results
171-
// ---------------------------------------------------------------------------
172-
PetscCall( PrintOutput(ceed, mem_type_backend,
173-
snes, ksp, U, l2_error_u, l2_error_p, app_ctx) );
176+
// ---------------------------------------------------------------------------
177+
// Print output results
178+
// ---------------------------------------------------------------------------
179+
PetscCall( PrintOutput(ceed, mem_type_backend,
180+
snes, ksp, U, l2_error_u, l2_error_p, app_ctx) );
174181

175-
// ---------------------------------------------------------------------------
176-
// Save solution (paraview)
177-
// ---------------------------------------------------------------------------
178-
PetscViewer viewer;
182+
// ---------------------------------------------------------------------------
183+
// Save solution (paraview)
184+
// ---------------------------------------------------------------------------
185+
PetscViewer viewer;
179186

180-
PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) );
181-
PetscCall( VecView(U, viewer) );
182-
PetscCall( PetscViewerDestroy(&viewer) );
187+
PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) );
188+
PetscCall( VecView(U, viewer) );
189+
PetscCall( PetscViewerDestroy(&viewer) );
183190

191+
PetscCall( SNESDestroy(&snes) );
192+
193+
}
184194
// ---------------------------------------------------------------------------
185195
// Free objects
186196
// ---------------------------------------------------------------------------
187197

188198
// Free PETSc objects
189199
PetscCall( DMDestroy(&dm) );
200+
PetscCall( DMDestroy(&dm_u0) );
201+
PetscCall( DMDestroy(&dm_p0) );
190202
PetscCall( VecDestroy(&U) );
191-
PetscCall( VecDestroy(&U_loc) );
192-
PetscCall( SNESDestroy(&snes) );
203+
PetscCall( CeedDataDestroy(ceed_data, problem_data) );
193204

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

197208
// -- Structs
198209
PetscCall( PetscFree(app_ctx) );
199210
PetscCall( PetscFree(problem_data) );
200-
PetscCall( PetscFree(phys_ctx) );
211+
PetscCall( PetscFree(ctx_initial_u0) );
212+
PetscCall( PetscFree(ctx_initial_p0) );
213+
PetscCall( PetscFree(ctx_residual_ut) );
201214

202215
// Free libCEED objects
203216
//CeedVectorDestroy(&bc_pressure);
204-
PetscCall( CeedDataDestroy(ceed_data) );
205217
CeedDestroy(&ceed);
206218

207219
return PetscFinalize();

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

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#include "../qfunctions/darcy-true2d.h"
2222
#include "../qfunctions/darcy-system2d.h"
2323
#include "../qfunctions/darcy-error2d.h"
24-
#include "../qfunctions/pressure-boundary2d.h"
24+
//#include "../qfunctions/pressure-boundary2d.h"
2525

2626
PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
2727
AppCtx app_ctx = *(AppCtx *)ctx;
@@ -47,8 +47,8 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
4747
problem_data->jacobian_loc = JacobianDarcySystem2D_loc;
4848
problem_data->error = DarcyError2D;
4949
problem_data->error_loc = DarcyError2D_loc;
50-
problem_data->bc_pressure = BCPressure2D;
51-
problem_data->bc_pressure_loc = BCPressure2D_loc;
50+
//problem_data->bc_pressure = BCPressure2D;
51+
//problem_data->bc_pressure_loc = BCPressure2D_loc;
5252
problem_data->has_ts = PETSC_FALSE;
5353

5454
// ------------------------------------------------------
@@ -77,9 +77,9 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
7777
CeedQFunctionContextCreate(ceed, &darcy_context);
7878
CeedQFunctionContextSetData(darcy_context, CEED_MEM_HOST, CEED_COPY_VALUES,
7979
sizeof(*darcy_ctx), darcy_ctx);
80-
problem_data->qfunction_context = darcy_context;
8180
CeedQFunctionContextSetDataDestroy(darcy_context, CEED_MEM_HOST,
8281
FreeContextPetsc);
83-
PetscCall( PetscFree(darcy_ctx) );
82+
problem_data->qfunction_context = darcy_context;
83+
8484
PetscFunctionReturn(0);
8585
}

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

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -21,7 +21,7 @@
2121
#include "../qfunctions/darcy-true3d.h"
2222
#include "../qfunctions/darcy-system3d.h"
2323
#include "../qfunctions/darcy-error3d.h"
24-
#include "../qfunctions/pressure-boundary3d.h"
24+
//#include "../qfunctions/pressure-boundary3d.h"
2525

2626
PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx) {
2727
AppCtx app_ctx = *(AppCtx *)ctx;
@@ -47,8 +47,8 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx) {
4747
problem_data->jacobian_loc = JacobianDarcySystem3D_loc;
4848
problem_data->error = DarcyError3D;
4949
problem_data->error_loc = DarcyError3D_loc;
50-
problem_data->bc_pressure = BCPressure3D;
51-
problem_data->bc_pressure_loc = BCPressure3D_loc;
50+
//problem_data->bc_pressure = BCPressure3D;
51+
//problem_data->bc_pressure_loc = BCPressure3D_loc;
5252
problem_data->has_ts = PETSC_FALSE;
5353

5454
// ------------------------------------------------------
@@ -77,9 +77,9 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx) {
7777
CeedQFunctionContextCreate(ceed, &darcy_context);
7878
CeedQFunctionContextSetData(darcy_context, CEED_MEM_HOST, CEED_COPY_VALUES,
7979
sizeof(*darcy_ctx), darcy_ctx);
80-
problem_data->qfunction_context = darcy_context;
8180
CeedQFunctionContextSetDataDestroy(darcy_context, CEED_MEM_HOST,
8281
FreeContextPetsc);
83-
PetscCall( PetscFree(darcy_ctx) );
82+
problem_data->qfunction_context = darcy_context;
83+
8484
PetscFunctionReturn(0);
8585
}

0 commit comments

Comments
 (0)