Skip to content

Commit 662499d

Browse files
rezgarshakerirezgarshakeri
committed
Added 3D version of Richards problem
1 parent 8ad7452 commit 662499d

File tree

8 files changed

+848
-0
lines changed

8 files changed

+848
-0
lines changed

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

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,8 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx);
1919
// 4) richard
2020
PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx);
2121

22+
PetscErrorCode Hdiv_RICHARD3D(Ceed ceed, ProblemData problem_data, void *ctx);
23+
2224
extern int FreeContextPetsc(void *);
2325

2426
#endif // register_problems_h

‎examples/Hdiv-mixed/main.c

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,8 @@
2626
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
2727
// ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4
2828
// ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
29+
// ./main -pc_type svd -problem richard2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
30+
// (boundary is not working)
2931
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
3032
// ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
3133
const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";

‎examples/Hdiv-mixed/problems/register-problem.c

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,8 @@ PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx) {
3434
// 4) richard
3535
PetscCall( PetscFunctionListAdd(&app_ctx->problems, "richard2d",
3636
Hdiv_RICHARD2D) );
37+
PetscCall( PetscFunctionListAdd(&app_ctx->problems, "richard3d",
38+
Hdiv_RICHARD3D) );
3739
PetscFunctionReturn(0);
3840
}
3941

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

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2+
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3+
// reserved. See files LICENSE and NOTICE for details.
4+
//
5+
// This file is part of CEED, a collection of benchmarks, miniapps, software
6+
// libraries and APIs for efficient high-order finite element and spectral
7+
// element discretizations for exascale applications. For more information and
8+
// source code availability see http://github.com/ceed.
9+
//
10+
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11+
// a collaborative effort of two U.S. Department of Energy organizations (Office
12+
// of Science and the National Nuclear Security Administration) responsible for
13+
// the planning and preparation of a capable exascale ecosystem, including
14+
// software, applications, hardware, advanced system engineering and early
15+
// testbed platforms, in support of the nation's exascale computing imperative.
16+
17+
/// @file
18+
/// Utility functions for setting up Richard problem in 3D
19+
20+
#include "../include/register-problem.h"
21+
#include "../qfunctions/richard-system3d.h"
22+
#include "../qfunctions/richard-true3d.h"
23+
#include "../qfunctions/richard-ics3d.h"
24+
#include "../qfunctions/darcy-error3d.h"
25+
#include "../qfunctions/post-processing3d.h"
26+
//#include "../qfunctions/pressure-boundary2d.h"
27+
#include "petscsystypes.h"
28+
29+
PetscErrorCode Hdiv_RICHARD3D(Ceed ceed, ProblemData problem_data, void *ctx) {
30+
AppCtx app_ctx = *(AppCtx *)ctx;
31+
RICHARDContext richard_ctx;
32+
CeedQFunctionContext richard_context;
33+
34+
PetscFunctionBeginUser;
35+
36+
PetscCall( PetscCalloc1(1, &richard_ctx) );
37+
38+
// ------------------------------------------------------
39+
// SET UP POISSON_QUAD2D
40+
// ------------------------------------------------------
41+
problem_data->dim = 3;
42+
problem_data->elem_node = 8;
43+
problem_data->q_data_size_face = 4;
44+
problem_data->quadrature_mode = CEED_GAUSS;
45+
problem_data->true_solution = RichardTrue3D;
46+
problem_data->true_solution_loc = RichardTrue3D_loc;
47+
problem_data->rhs_u0 = RichardRhsU03D;
48+
problem_data->rhs_u0_loc = RichardRhsU03D_loc;
49+
problem_data->ics_u = RichardICsU3D;
50+
problem_data->ics_u_loc = RichardICsU3D_loc;
51+
problem_data->rhs_p0 = RichardRhsP03D;
52+
problem_data->rhs_p0_loc = RichardRhsP03D_loc;
53+
problem_data->ics_p = RichardICsP3D;
54+
problem_data->ics_p_loc = RichardICsP3D_loc;
55+
problem_data->residual = RichardSystem3D;
56+
problem_data->residual_loc = RichardSystem3D_loc;
57+
//problem_data->jacobian = JacobianRichardSystem2D;
58+
//problem_data->jacobian_loc = JacobianRichardSystem2D_loc;
59+
problem_data->error = DarcyError3D;
60+
problem_data->error_loc = DarcyError3D_loc;
61+
//problem_data->bc_pressure = BCPressure2D;
62+
//problem_data->bc_pressure_loc = BCPressure2D_loc;
63+
problem_data->post_rhs = PostProcessingRhs3D;
64+
problem_data->post_rhs_loc = PostProcessingRhs3D_loc;
65+
problem_data->post_mass = PostProcessingMass3D;
66+
problem_data->post_mass_loc = PostProcessingMass3D_loc;
67+
problem_data->has_ts = PETSC_TRUE;
68+
problem_data->view_solution = app_ctx->view_solution;
69+
70+
// ------------------------------------------------------
71+
// Command line Options
72+
// ------------------------------------------------------
73+
CeedScalar kappa = 10., alpha_a = 1., b_a = 10., rho_a0 = 998.2,
74+
beta = 0., g = 9.8, p0 = 101325;
75+
76+
PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
77+
PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
78+
kappa, &kappa, NULL));
79+
PetscCall( PetscOptionsScalar("-alpha_a", "Parameter for relative permeability",
80+
NULL,
81+
alpha_a, &alpha_a, NULL));
82+
PetscCall( PetscOptionsScalar("-b_a", "Parameter for relative permeability",
83+
NULL,
84+
b_a, &b_a, NULL));
85+
PetscCall( PetscOptionsScalar("-rho_a0", "Density at p0", NULL,
86+
rho_a0, &rho_a0, NULL));
87+
PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
88+
beta, &beta, NULL));
89+
app_ctx->t_final = 0.5;
90+
PetscCall( PetscOptionsScalar("-t_final", "End time", NULL,
91+
app_ctx->t_final, &app_ctx->t_final, NULL));
92+
PetscOptionsEnd();
93+
94+
richard_ctx->kappa = kappa;
95+
richard_ctx->alpha_a = alpha_a;
96+
richard_ctx->b_a = b_a;
97+
richard_ctx->rho_a0 = rho_a0;
98+
richard_ctx->beta = beta;
99+
richard_ctx->g = g;
100+
richard_ctx->p0 = p0;
101+
richard_ctx->gamma = 5.;
102+
richard_ctx->t = 0.;
103+
richard_ctx->t_final = app_ctx->t_final;
104+
CeedQFunctionContextCreate(ceed, &richard_context);
105+
CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
106+
sizeof(*richard_ctx), richard_ctx);
107+
//CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
108+
// FreeContextPetsc);
109+
CeedQFunctionContextRegisterDouble(richard_context, "time",
110+
offsetof(struct RICHARDContext_, t), 1, "current solver time");
111+
CeedQFunctionContextRegisterDouble(richard_context, "final_time",
112+
offsetof(struct RICHARDContext_, t_final), 1, "final time");
113+
CeedQFunctionContextRegisterDouble(richard_context, "time_step",
114+
offsetof(struct RICHARDContext_, dt), 1, "time step");
115+
problem_data->true_qfunction_ctx = richard_context;
116+
CeedQFunctionContextReferenceCopy(richard_context,
117+
&problem_data->rhs_u0_qfunction_ctx);
118+
CeedQFunctionContextReferenceCopy(richard_context,
119+
&problem_data->residual_qfunction_ctx);
120+
CeedQFunctionContextReferenceCopy(richard_context,
121+
&problem_data->error_qfunction_ctx);
122+
PetscCall( PetscFree(richard_ctx) );
123+
124+
PetscFunctionReturn(0);
125+
}

‎examples/Hdiv-mixed/qfunctions/post-processing3d.h

Lines changed: 100 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,100 @@
1+
// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2+
// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3+
// reserved. See files LICENSE and NOTICE for details.
4+
//
5+
// This file is part of CEED, a collection of benchmarks, miniapps, software
6+
// libraries and APIs for efficient high-order finite element and spectral
7+
// element discretizations for exascale applications. For more information and
8+
// source code availability see http://github.com/ceed.
9+
//
10+
// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
11+
// a collaborative effort of two U.S. Department of Energy organizations (Office
12+
// of Science and the National Nuclear Security Administration) responsible for
13+
// the planning and preparation of a capable exascale ecosystem, including
14+
// software, applications, hardware, advanced system engineering and early
15+
// testbed platforms, in support of the nation's exascale computing imperative.
16+
17+
/// @file
18+
/// Force of Richard problem 2D (quad element) using PETSc
19+
20+
#ifndef POST_PROCESSING3D_H
21+
#define POST_PROCESSING3D_H
22+
23+
#include <math.h>
24+
#include <ceed.h>
25+
#include "ceed/ceed-f64.h"
26+
#include "utils.h"
27+
28+
// -----------------------------------------------------------------------------
29+
// We solve (v, u) = (v, uh), to project Hdiv to L2 space
30+
// This QFunction create post_rhs = (v, uh)
31+
// -----------------------------------------------------------------------------
32+
CEED_QFUNCTION(PostProcessingRhs3D)(void *ctx, const CeedInt Q,
33+
const CeedScalar *const *in,
34+
CeedScalar *const *out) {
35+
// *INDENT-OFF*
36+
// Inputs
37+
const CeedScalar (*w) = in[0],
38+
(*dxdX)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1],
39+
(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
40+
// Outputs
41+
CeedScalar (*post_rhs) = out[0];
42+
43+
// Quadrature Point Loop
44+
CeedPragmaSIMD
45+
for (CeedInt i=0; i<Q; i++) {
46+
const CeedScalar J[3][3] = {{dxdX[0][0][i], dxdX[1][0][i], dxdX[2][0][i]},
47+
{dxdX[0][1][i], dxdX[1][1][i], dxdX[2][1][i]},
48+
{dxdX[0][2][i], dxdX[1][2][i], dxdX[2][2][i]}};
49+
50+
// 1) Compute Piola map: uh = J*u/detJ
51+
// 2) rhs = (v, uh) = uh*w*det_J
52+
// ==> rhs = J*u*w
53+
CeedScalar u1[3] = {u[0][i], u[1][i], u[2][i]}, rhs[3];
54+
AlphaMatVecMult3x3(w[i], J, u1, rhs);
55+
56+
post_rhs[i+0*Q] = rhs[0];
57+
post_rhs[i+1*Q] = rhs[1];
58+
post_rhs[i+2*Q] = rhs[2];
59+
} // End of Quadrature Point Loop
60+
return 0;
61+
}
62+
63+
// -----------------------------------------------------------------------------
64+
// We solve (v, u) = (v, uh), to project Hdiv to L2 space
65+
// This QFunction create mass matrix (v, u), then we solve using ksp to have
66+
// projected uh in L2 space and use it for post-processing
67+
// -----------------------------------------------------------------------------
68+
CEED_QFUNCTION(PostProcessingMass3D)(void *ctx, const CeedInt Q,
69+
const CeedScalar *const *in,
70+
CeedScalar *const *out) {
71+
// *INDENT-OFF*
72+
// Inputs
73+
const CeedScalar (*w) = in[0],
74+
(*dxdX)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1],
75+
(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2];
76+
77+
// Outputs
78+
CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0];
79+
80+
// Quadrature Point Loop
81+
CeedPragmaSIMD
82+
for (CeedInt i=0; i<Q; i++) {
83+
// *INDENT-OFF*
84+
// Setup, J = dx/dX
85+
const CeedScalar J[3][3] = {{dxdX[0][0][i], dxdX[1][0][i], dxdX[2][0][i]},
86+
{dxdX[0][1][i], dxdX[1][1][i], dxdX[2][1][i]},
87+
{dxdX[0][2][i], dxdX[1][2][i], dxdX[2][2][i]}};
88+
const CeedScalar det_J = MatDet3x3(J);
89+
90+
// *INDENT-ON*
91+
// (v, u): v = u*w*detJ
92+
for (CeedInt k = 0; k < 3; k++) {
93+
v[k][i] = u[k][i]*w[i]*det_J;
94+
}
95+
} // End of Quadrature Point Loop
96+
return 0;
97+
}
98+
// -----------------------------------------------------------------------------
99+
100+
#endif //End of POST_PROCESSING3D_H

0 commit comments

Comments
 (0)