Skip to content

Commit 4ebe4aa

Browse files
rezgarshakerirezgarshakeri
committed
Create Initial conditions for Richards problem
1 parent 833e410 commit 4ebe4aa

12 files changed

Lines changed: 411 additions & 332 deletions

File tree

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

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,8 @@
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 U_loc,
10+
Vec U);
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,

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

Lines changed: 7 additions & 46 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,6 +31,8 @@ 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
@@ -81,10 +41,11 @@ 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,
8343
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;
44+
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics, qf_true;
45+
CeedOperator op_residual, op_jacobian, op_error, op_ics, op_true;
46+
CeedVector x_ceed, y_ceed, x_coord, x0_ceed, x_t_ceed;
8747
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut;
48+
CeedInt num_elem;
8849
};
8950

9051
// Problem specific data

‎examples/Hdiv-mixed/main.c‎

Lines changed: 46 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -73,9 +73,9 @@ 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;
77+
PetscCall( PetscCalloc1(1, &ctx_residual_ut) );
78+
ceed_data->ctx_residual_ut = ctx_residual_ut;
7979
// ---------------------------------------------------------------------------
8080
// Process command line options
8181
// ---------------------------------------------------------------------------
@@ -106,7 +106,7 @@ int main(int argc, char **argv) {
106106
PetscCall( CreateDM(comm, vec_type, &dm) );
107107
// TODO: add mesh option
108108
// perturb to have smooth random mesh
109-
// PetscCall( PerturbVerticesSmooth(dm) );
109+
PetscCall( PerturbVerticesSmooth(dm) );
110110

111111
// ---------------------------------------------------------------------------
112112
// Setup FE
@@ -116,10 +116,12 @@ int main(int argc, char **argv) {
116116
// ---------------------------------------------------------------------------
117117
// Create local Force vector
118118
// ---------------------------------------------------------------------------
119-
Vec U_loc;
119+
Vec U_loc, U;
120120
PetscInt U_loc_size;
121-
//CeedVector bc_pressure;
121+
PetscCall( DMCreateGlobalVector(dm, &U) );
122+
PetscCall( VecZeroEntries(U) );
122123
PetscCall( DMCreateLocalVector(dm, &U_loc) );
124+
PetscCall( VecZeroEntries(U_loc) );
123125
// Local size for libCEED
124126
PetscCall( VecGetSize(U_loc, &U_loc_size) );
125127

@@ -141,46 +143,47 @@ int main(int argc, char **argv) {
141143
// ---------------------------------------------------------------------------
142144
// Create global initial conditions
143145
// ---------------------------------------------------------------------------
144-
Vec U0;
145-
CreateInitialConditions(dm, ceed_data, &U0);
146-
VecView(U0, PETSC_VIEWER_STDOUT_WORLD);
147-
PetscCall( VecDestroy(&U0) );
146+
CreateInitialConditions(dm, ceed_data, U_loc, U);
147+
VecView(U, PETSC_VIEWER_STDOUT_WORLD);
148148
}
149149

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);
150+
if (!problem_data->has_ts) {
151+
// ---------------------------------------------------------------------------
152+
// Solve PDE
153+
// ---------------------------------------------------------------------------
154+
// Create SNES
155+
SNES snes;
156+
KSP ksp;
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);
161161

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) );
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) );
168168

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) );
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) );
174174

175-
// ---------------------------------------------------------------------------
176-
// Save solution (paraview)
177-
// ---------------------------------------------------------------------------
178-
PetscViewer viewer;
175+
// ---------------------------------------------------------------------------
176+
// Save solution (paraview)
177+
// ---------------------------------------------------------------------------
178+
PetscViewer viewer;
179179

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

184+
PetscCall( SNESDestroy(&snes) );
185+
186+
}
184187
// ---------------------------------------------------------------------------
185188
// Free objects
186189
// ---------------------------------------------------------------------------
@@ -189,19 +192,20 @@ int main(int argc, char **argv) {
189192
PetscCall( DMDestroy(&dm) );
190193
PetscCall( VecDestroy(&U) );
191194
PetscCall( VecDestroy(&U_loc) );
192-
PetscCall( SNESDestroy(&snes) );
193195

194196
// -- Function list
195197
PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );
196198

197199
// -- Structs
198200
PetscCall( PetscFree(app_ctx) );
199201
PetscCall( PetscFree(problem_data) );
200-
PetscCall( PetscFree(phys_ctx) );
202+
PetscCall( PetscFree(ctx_residual_ut) );
201203

202204
// Free libCEED objects
203205
//CeedVectorDestroy(&bc_pressure);
204-
PetscCall( CeedDataDestroy(ceed_data) );
206+
//PetscCall( CeedDataDestroy(ceed_data) );
207+
CeedQFunctionDestroy(&ceed_data->qf_true);
208+
CeedOperatorDestroy(&ceed_data->op_true);
205209
CeedDestroy(&ceed);
206210

207211
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
}

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

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

2626
PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
@@ -49,8 +49,8 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
4949
//problem_data->jacobian_loc = JacobianRichardSystem2D_loc;
5050
//problem_data->error = DarcyError2D;
5151
//problem_data->error_loc = DarcyError2D_loc;
52-
problem_data->bc_pressure = BCPressure2D;
53-
problem_data->bc_pressure_loc = BCPressure2D_loc;
52+
//problem_data->bc_pressure = BCPressure2D;
53+
//problem_data->bc_pressure_loc = BCPressure2D_loc;
5454
problem_data->has_ts = PETSC_TRUE;
5555

5656
// ------------------------------------------------------
@@ -72,6 +72,9 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
7272
rho_a0, &rho_a0, NULL));
7373
PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
7474
beta, &beta, NULL));
75+
app_ctx->t_final = 5.;
76+
PetscCall( PetscOptionsScalar("-t_final", "End time", NULL,
77+
app_ctx->t_final, &app_ctx->t_final, NULL));
7578
PetscOptionsEnd();
7679

7780
richard_ctx->kappa = kappa;
@@ -82,15 +85,18 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
8285
richard_ctx->g = g;
8386
richard_ctx->p0 = p0;
8487
richard_ctx->gamma = 5.;
88+
richard_ctx->t = 0.;
89+
richard_ctx->t_final = app_ctx->t_final;
8590
CeedQFunctionContextCreate(ceed, &richard_context);
8691
CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
8792
sizeof(*richard_ctx), richard_ctx);
93+
CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
94+
FreeContextPetsc);
8895
CeedQFunctionContextRegisterDouble(richard_context, "time",
8996
offsetof(struct RICHARDContext_, t), 1, "current solver time");
97+
CeedQFunctionContextRegisterDouble(richard_context, "final_time",
98+
offsetof(struct RICHARDContext_, t_final), 1, "final time");
9099
problem_data->qfunction_context = richard_context;
91-
CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
92-
FreeContextPetsc);
93100

94-
PetscCall( PetscFree(richard_ctx) );
95101
PetscFunctionReturn(0);
96102
}

‎examples/Hdiv-mixed/qfunctions/richard-system2d.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -71,7 +71,7 @@ struct RICHARDContext_ {
7171
CeedScalar rho_a0;
7272
CeedScalar alpha_a, b_a;
7373
CeedScalar beta, p0;
74-
CeedScalar t;
74+
CeedScalar t, t_final;
7575
CeedScalar gamma;
7676
};
7777
#endif

0 commit comments

Comments
 (0)