Skip to content

Commit bfad0b5

Browse files
rezgarshakerirezgarshakeri
committed
setup qfunction context for darcy problem
1 parent 3569188 commit bfad0b5

23 files changed

Lines changed: 991 additions & 405 deletions

‎examples/Hdiv-mixed/include/register-problem.h‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -9,13 +9,15 @@ PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx);
99
// Set up problems function prototype
1010
// -----------------------------------------------------------------------------
1111
// 1) darcy2d
12-
PetscErrorCode Hdiv_DARCY2D(ProblemData problem_data, void *ctx);
12+
PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx);
1313

1414
// 2) darcy3d
15-
PetscErrorCode Hdiv_DARCY3D(ProblemData problem_data, void *ctx);
15+
PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx);
1616

1717
// 3) darcy3dprism
1818

1919
// 4) richard
2020

21+
extern int FreeContextPetsc(void *);
22+
2123
#endif // register_problems_h

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

Lines changed: 2 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8,9 +8,8 @@
88
#include "structs.h"
99

1010
// ---------------------------------------------------------------------------
11-
// Set-up DM
11+
// Setup DM
1212
// ---------------------------------------------------------------------------
13-
PetscErrorCode CreateDistributedDM(MPI_Comm comm, ProblemData problem_data,
14-
DM *dm);
13+
PetscErrorCode CreateDM(MPI_Comm comm, VecType vec_type, DM *dm);
1514

1615
#endif // setupdm_h

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

Lines changed: 15 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,15 @@
1+
#ifndef setupfe_h
2+
#define setupfe_h
3+
4+
#include <petsc.h>
5+
#include <petscdmplex.h>
6+
#include <petscsys.h>
7+
#include <ceed.h>
8+
#include "structs.h"
9+
10+
// ---------------------------------------------------------------------------
11+
// Setup FE
12+
// ---------------------------------------------------------------------------
13+
PetscErrorCode SetupFE(MPI_Comm comm, DM dm);
14+
15+
#endif // setupfe_h

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

Lines changed: 12 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -6,24 +6,21 @@
66

77
#include "structs.h"
88

9-
PetscErrorCode SetupCommonCtx(DM dm, Ceed ceed,
10-
CeedData ceed_data,
11-
OperatorApplyContext op_apply_ctx);
12-
PetscErrorCode SetupJacobianOperatorCtx(CeedData ceed_data,
13-
OperatorApplyContext op_apply_ctx);
14-
PetscErrorCode SetupResidualOperatorCtx(CeedData ceed_data,
15-
OperatorApplyContext op_apply_ctx);
16-
PetscErrorCode SetupMMSOperatorCtx(CeedData ceed_data,
17-
OperatorApplyContext op_apply_ctx);
9+
PetscErrorCode SetupJacobianOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
10+
OperatorApplyContext ctx_jacobian);
11+
PetscErrorCode SetupResidualOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
12+
OperatorApplyContext ctx_residual);
13+
PetscErrorCode SetupErrorOperatorCtx(DM dm, Ceed ceed, CeedData ceed_data,
14+
OperatorApplyContext ctx_error);
1815
PetscErrorCode ApplyJacobian(Mat A, Vec X, Vec Y);
1916
PetscErrorCode SNESFormResidual(SNES snes, Vec X, Vec Y, void *ctx);
2017
PetscErrorCode SNESFormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx);
21-
PetscErrorCode PDESolver(CeedData ceed_data, VecType vec_type, SNES snes,
22-
KSP ksp,
23-
Vec F, Vec *U_g, OperatorApplyContext op_apply_ctx);
24-
PetscErrorCode ComputeL2Error(CeedData ceed_data, Vec U, CeedVector target,
25-
CeedScalar *l2_error_u, CeedScalar *l2_error_p,
26-
OperatorApplyContext op_apply_ctx);
18+
PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data,
19+
VecType vec_type, SNES snes, KSP ksp,
20+
Vec F, Vec *U_g);
21+
PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U,
22+
CeedVector target,
23+
CeedScalar *l2_error_u, CeedScalar *l2_error_p);
2724
PetscErrorCode PrintOutput(Ceed ceed,
2825
CeedMemType mem_type_backend,
2926
SNES snes, KSP ksp,

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

Lines changed: 44 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -18,63 +18,49 @@ struct AppCtx_ {
1818
// Problem type arguments
1919
PetscFunctionList problems;
2020
char problem_name[PETSC_MAX_PATH_LEN];
21-
2221
};
2322

24-
// libCEED data struct
25-
typedef struct CeedData_ *CeedData;
26-
struct CeedData_ {
27-
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
28-
CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_U_i,
29-
elem_restr_p;
30-
CeedQFunction qf_residual, qf_jacobian, qf_error;
31-
CeedOperator op_residual, op_jacobian, op_error;
32-
CeedVector x_ceed, y_ceed, x_coord;
33-
};
23+
// 2) richard
24+
// We have 3 experiment parameters as described in Table 1:P1, P2, P3
25+
// Matthew Farthing, Christopher Kees, Cass Miller (2003)
26+
// https://www.sciencedirect.com/science/article/pii/S0309170802001872
3427

35-
// 1) darcy2d
36-
#ifndef PHYSICS_DARCY2D_STRUCT
37-
#define PHYSICS_DARCY2D_STRUCT
38-
typedef struct DARCY2DContext_ *DARCY2DContext;
39-
struct DARCY2DContext_ {
40-
CeedScalar kappa;
28+
#ifndef PHYSICS_RICHARDP2_STRUCT
29+
#define PHYSICS_RICHARDP2_STRUCT
30+
typedef struct RICHARDP2Context_ *RICHARDP2Context;
31+
struct RICHARDP2Context_ {
32+
CeedScalar K_star;
33+
CeedScalar theta_s;
34+
CeedScalar theta_r;
35+
CeedScalar alpha_v;
36+
CeedScalar n_v;
37+
CeedScalar m_v;
38+
CeedScalar m_r;
39+
CeedScalar rho_0;
40+
CeedScalar beta;
4141
};
4242
#endif
4343

44-
// 2) darcy3d
45-
#ifndef PHYSICS_DARCY3D_STRUCT
46-
#define PHYSICS_DARCY3D_STRUCT
47-
typedef struct DARCY3DContext_ *DARCY3DContext;
48-
struct DARCY3DContext_ {
49-
CeedScalar kappa;
50-
};
51-
#endif
52-
53-
// 3) richard2d
54-
#ifndef PHYSICS_RICHARD2D_STRUCT
55-
#define PHYSICS_RICHARD2D_STRUCT
56-
typedef struct RICHARD2DContext_ *RICHARD2DContext;
57-
struct RICHARD2DContext_ {
58-
CeedScalar kappa;
59-
};
60-
#endif
61-
62-
// 4) richard3d
63-
#ifndef PHYSICS_RICHARD3D_STRUCT
64-
#define PHYSICS_RICHARD3D_STRUCT
65-
typedef struct RICHARD3DContext_ *RICHARD3DContext;
66-
struct RICHARD3DContext_ {
67-
CeedScalar kappa;
44+
#ifndef PHYSICS_RICHARDP3_STRUCT
45+
#define PHYSICS_RICHARDP3_STRUCT
46+
typedef struct RICHARDP3Context_ *RICHARDP3Context;
47+
struct RICHARDP3Context_ {
48+
CeedScalar K_star;
49+
CeedScalar theta_s;
50+
CeedScalar theta_r;
51+
CeedScalar alpha_star_v;
52+
CeedScalar n_v;
53+
CeedScalar m_v;
54+
CeedScalar rho_0;
55+
CeedScalar beta;
6856
};
6957
#endif
7058

7159
// Struct that contains all enums and structs used for the physics of all problems
7260
typedef struct Physics_ *Physics;
7361
struct Physics_ {
74-
DARCY2DContext darcy2d_ctx;
75-
DARCY3DContext darcy3d_ctx;
76-
RICHARD2DContext richard2d_ctx;
77-
RICHARD3DContext richard3d_ctx;
62+
RICHARDP2Context richard_p2_ctx;
63+
RICHARDP3Context richard_p3_ctx;
7864
};
7965

8066
// PETSc operator contexts
@@ -88,6 +74,18 @@ struct OperatorApplyContext_ {
8874
Ceed ceed;
8975
};
9076

77+
// libCEED data struct
78+
typedef struct CeedData_ *CeedData;
79+
struct CeedData_ {
80+
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
81+
CeedElemRestriction elem_restr_x, elem_restr_u, elem_restr_U_i,
82+
elem_restr_p;
83+
CeedQFunction qf_residual, qf_jacobian, qf_error;
84+
CeedOperator op_residual, op_jacobian, op_error;
85+
CeedVector x_ceed, y_ceed, x_coord;
86+
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error;
87+
};
88+
9189
// Problem specific data
9290
typedef struct ProblemData_ *ProblemData;
9391
struct ProblemData_ {
@@ -97,7 +95,7 @@ struct ProblemData_ {
9795
*error_loc, *setup_true_loc, *bc_pressure_loc;
9896
CeedQuadMode quadrature_mode;
9997
CeedInt elem_node, dim, q_data_size_face;
100-
98+
CeedQFunctionContext qfunction_context;
10199
};
102100

103101
#endif // structs_h

‎examples/Hdiv-mixed/main.c‎

Lines changed: 34 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -38,6 +38,29 @@ int main(int argc, char **argv) {
3838
// ---------------------------------------------------------------------------
3939
PetscCall( PetscInitialize(&argc, &argv, NULL, help) );
4040

41+
// ---------------------------------------------------------------------------
42+
// Initialize libCEED
43+
// ---------------------------------------------------------------------------
44+
// -- Initialize backend
45+
Ceed ceed;
46+
CeedInit("/cpu/self/ref/serial", &ceed);
47+
CeedMemType mem_type_backend;
48+
CeedGetPreferredMemType(ceed, &mem_type_backend);
49+
50+
VecType vec_type = NULL;
51+
switch (mem_type_backend) {
52+
case CEED_MEM_HOST: vec_type = VECSTANDARD; break;
53+
case CEED_MEM_DEVICE: {
54+
const char *resolved;
55+
CeedGetResource(ceed, &resolved);
56+
if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
57+
else if (strstr(resolved, "/gpu/hip/occa"))
58+
vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678
59+
else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
60+
else vec_type = VECSTANDARD;
61+
}
62+
}
63+
4164
// ---------------------------------------------------------------------------
4265
// Create structs
4366
// ---------------------------------------------------------------------------
@@ -71,43 +94,24 @@ int main(int argc, char **argv) {
7194
// Choose the problem from the list of registered problems
7295
// ---------------------------------------------------------------------------
7396
{
74-
PetscErrorCode (*p)(ProblemData, void *);
97+
PetscErrorCode (*p)(Ceed, ProblemData, void *);
7598
PetscCall( PetscFunctionListFind(app_ctx->problems, app_ctx->problem_name,
7699
&p) );
77100
if (!p) SETERRQ(PETSC_COMM_SELF, 1, "Problem '%s' not found",
78101
app_ctx->problem_name);
79-
PetscCall( (*p)(problem_data, &phys_ctx) );
102+
PetscCall( (*p)(ceed, problem_data, &app_ctx) );
80103
}
81104

82105
// ---------------------------------------------------------------------------
83-
// Initialize libCEED
84-
// ---------------------------------------------------------------------------
85-
// -- Initialize backend
86-
Ceed ceed;
87-
CeedInit("/cpu/self/ref/serial", &ceed);
88-
CeedMemType mem_type_backend;
89-
CeedGetPreferredMemType(ceed, &mem_type_backend);
90-
// ---------------------------------------------------------------------------
91-
// Setup DM
106+
// Create DM
92107
// ---------------------------------------------------------------------------
93-
// PETSc objects
94108
DM dm;
95-
VecType vec_type = NULL;
96-
PetscCall( CreateDistributedDM(comm, problem_data, &dm) );
97-
switch (mem_type_backend) {
98-
case CEED_MEM_HOST: vec_type = VECSTANDARD; break;
99-
case CEED_MEM_DEVICE: {
100-
const char *resolved;
101-
CeedGetResource(ceed, &resolved);
102-
if (strstr(resolved, "/gpu/cuda")) vec_type = VECCUDA;
103-
else if (strstr(resolved, "/gpu/hip/occa"))
104-
vec_type = VECSTANDARD; // https://github.com/CEED/libCEED/issues/678
105-
else if (strstr(resolved, "/gpu/hip")) vec_type = VECHIP;
106-
else vec_type = VECSTANDARD;
107-
}
108-
}
109-
PetscCall( DMSetVecType(dm, vec_type) );
109+
PetscCall( CreateDM(comm, vec_type, &dm) );
110110

111+
// ---------------------------------------------------------------------------
112+
// Setup FE
113+
// ---------------------------------------------------------------------------
114+
SetupFE(comm, dm);
111115

112116
// ---------------------------------------------------------------------------
113117
// Create local Force vector
@@ -153,19 +157,17 @@ int main(int argc, char **argv) {
153157
SNES snes;
154158
KSP ksp;
155159
Vec U;
156-
op_apply_ctx->comm = comm;
157160
PetscCall( SNESCreate(comm, &snes) );
158161
PetscCall( SNESGetKSP(snes, &ksp) );
159-
PetscCall( SetupCommonCtx(dm, ceed, ceed_data, op_apply_ctx) );
160-
PetscCall( PDESolver(ceed_data, vec_type, snes, ksp, F, &U, op_apply_ctx) );
162+
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, F, &U) );
161163
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
162164

163165
// ---------------------------------------------------------------------------
164166
// Compute L2 error of mms problem
165167
// ---------------------------------------------------------------------------
166168
CeedScalar l2_error_u, l2_error_p;
167-
PetscCall( ComputeL2Error(ceed_data, U, target, &l2_error_u, &l2_error_p,
168-
op_apply_ctx) );
169+
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, target, &l2_error_u,
170+
&l2_error_p) );
169171

170172
// ---------------------------------------------------------------------------
171173
// Print output results
@@ -201,8 +203,6 @@ int main(int argc, char **argv) {
201203
// -- Structs
202204
PetscCall( PetscFree(app_ctx) );
203205
PetscCall( PetscFree(problem_data) );
204-
PetscCall( PetscFree(phys_ctx->darcy2d_ctx) );
205-
PetscCall( PetscFree(phys_ctx->darcy3d_ctx) );
206206
PetscCall( PetscFree(phys_ctx) );
207207
PetscCall( PetscFree(op_apply_ctx) );
208208

‎examples/Hdiv-mixed/main.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44

55
#include "include/setup-libceed.h"
66
#include "include/setup-dm.h"
7+
#include "include/setup-fe.h"
78
#include "include/setup-boundary.h"
89
#include "include/cl-options.h"
910
#include "include/register-problem.h"

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

Lines changed: 23 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -19,17 +19,18 @@
1919

2020
#include "../include/register-problem.h"
2121
#include "../qfunctions/darcy-force2d.h"
22-
#include "../qfunctions/darcy-mass2d.h"
22+
#include "../qfunctions/darcy-system2d.h"
2323
#include "../qfunctions/darcy-error2d.h"
2424
#include "../qfunctions/pressure-boundary2d.h"
2525

26-
// Hdiv_DARCY2D is registered in cl-option.c
27-
PetscErrorCode Hdiv_DARCY2D(ProblemData problem_data, void *ctx) {
28-
Physics phys = *(Physics *)ctx;
29-
MPI_Comm comm = PETSC_COMM_WORLD;
26+
PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
27+
AppCtx app_ctx = *(AppCtx *)ctx;
28+
DARCYContext darcy_ctx;
29+
CeedQFunctionContext darcy_context;
30+
3031
PetscFunctionBeginUser;
3132

32-
PetscCall( PetscCalloc1(1, &phys->darcy2d_ctx) );
33+
PetscCall( PetscCalloc1(1, &darcy_ctx) );
3334

3435
// ------------------------------------------------------
3536
// SET UP POISSON_QUAD2D
@@ -40,10 +41,10 @@ PetscErrorCode Hdiv_DARCY2D(ProblemData problem_data, void *ctx) {
4041
problem_data->quadrature_mode = CEED_GAUSS;
4142
problem_data->force = DarcyForce2D;
4243
problem_data->force_loc = DarcyForce2D_loc;
43-
problem_data->residual = DarcyMass2D;
44-
problem_data->residual_loc = DarcyMass2D_loc;
45-
problem_data->jacobian = JacobianDarcyMass2D;
46-
problem_data->jacobian_loc = JacobianDarcyMass2D_loc;
44+
problem_data->residual = DarcySystem2D;
45+
problem_data->residual_loc = DarcySystem2D_loc;
46+
problem_data->jacobian = JacobianDarcySystem2D;
47+
problem_data->jacobian_loc = JacobianDarcySystem2D_loc;
4748
problem_data->error = DarcyError2D;
4849
problem_data->error_loc = DarcyError2D_loc;
4950
problem_data->bc_pressure = BCPressure2D;
@@ -52,10 +53,19 @@ PetscErrorCode Hdiv_DARCY2D(ProblemData problem_data, void *ctx) {
5253
// ------------------------------------------------------
5354
// Command line Options
5455
// ------------------------------------------------------
55-
PetscOptionsBegin(comm, NULL, "Options for Hdiv-mixed problem", NULL);
56-
56+
CeedScalar kappa = 1.;
57+
PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
58+
PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
59+
kappa, &kappa, NULL));
5760
PetscOptionsEnd();
5861

59-
PetscCall( PetscFree(phys->darcy2d_ctx) );
62+
darcy_ctx->kappa = kappa;
63+
CeedQFunctionContextCreate(ceed, &darcy_context);
64+
CeedQFunctionContextSetData(darcy_context, CEED_MEM_HOST, CEED_COPY_VALUES,
65+
sizeof(*darcy_ctx), darcy_ctx);
66+
problem_data->qfunction_context = darcy_context;
67+
CeedQFunctionContextSetDataDestroy(darcy_context, CEED_MEM_HOST,
68+
FreeContextPetsc);
69+
//PetscCall( PetscFree(darcy_ctx) );
6070
PetscFunctionReturn(0);
6171
}

0 commit comments

Comments
 (0)