Skip to content

Commit 5e0f655

Browse files
committed
setup richards qfunction
1 parent 4b46c93 commit 5e0f655

File tree

7 files changed

+273
-79
lines changed

7 files changed

+273
-79
lines changed

examples/Hdiv-mixed/problems/richard2d.c

Lines changed: 35 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -19,16 +19,17 @@
1919

2020
#include "../include/register-problem.h"
2121
#include "../qfunctions/richard-system2d.h"
22+
#include "../qfunctions/richard-force2d.h"
2223
#include "../qfunctions/pressure-boundary2d.h"
2324

2425
PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
2526
AppCtx app_ctx = *(AppCtx *)ctx;
26-
//RICHARDContext richard_ctx;
27-
//CeedQFunctionContext richard_context;
27+
RICHARDContext richard_ctx;
28+
CeedQFunctionContext richard_context;
2829

2930
PetscFunctionBeginUser;
3031

31-
//PetscCall( PetscCalloc1(1, &richard_ctx) );
32+
PetscCall( PetscCalloc1(1, &richard_ctx) );
3233

3334
// ------------------------------------------------------
3435
// SET UP POISSON_QUAD2D
@@ -37,8 +38,8 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
3738
problem_data->elem_node = 4;
3839
problem_data->q_data_size_face = 3;
3940
problem_data->quadrature_mode = CEED_GAUSS;
40-
//problem_data->force = DarcyForce2D;
41-
//problem_data->force_loc = DarcyForce2D_loc;
41+
problem_data->force = RichardForce2D;
42+
problem_data->force_loc = RichardForce2D_loc;
4243
problem_data->residual = RichardSystem2D;
4344
problem_data->residual_loc = RichardSystem2D_loc;
4445
problem_data->jacobian = JacobianRichardSystem2D;
@@ -51,9 +52,38 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
5152
// ------------------------------------------------------
5253
// Command line Options
5354
// ------------------------------------------------------
55+
CeedScalar kappa = 10., alpha_a = 1., b_a = 10., rho_0 = 998.2,
56+
beta = 0., g = 9.8, p0 = 101325;
57+
5458
PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
59+
PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
60+
kappa, &kappa, NULL));
61+
PetscCall( PetscOptionsScalar("-alpha_a", "Parameter for relative permeability",
62+
NULL,
63+
alpha_a, &alpha_a, NULL));
64+
PetscCall( PetscOptionsScalar("-b_a", "Parameter for relative permeability",
65+
NULL,
66+
b_a, &b_a, NULL));
67+
PetscCall( PetscOptionsScalar("-rho_0", "Density at p0", NULL,
68+
rho_0, &rho_0, NULL));
69+
PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL,
70+
beta, &beta, NULL));
5571
PetscOptionsEnd();
5672

73+
richard_ctx->kappa = kappa;
74+
richard_ctx->alpha_a = alpha_a;
75+
richard_ctx->b_a = b_a;
76+
richard_ctx->rho_0 = rho_0;
77+
richard_ctx->beta = beta;
78+
richard_ctx->g = g;
79+
richard_ctx->p0 = p0;
80+
CeedQFunctionContextCreate(ceed, &richard_context);
81+
CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
82+
sizeof(*richard_ctx), richard_ctx);
83+
problem_data->qfunction_context = richard_context;
84+
CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
85+
FreeContextPetsc);
86+
5787
//PetscCall( PetscFree(richard_ctx) );
5888
PetscFunctionReturn(0);
5989
}

examples/Hdiv-mixed/qfunctions/darcy-system2d.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -184,4 +184,4 @@ CEED_QFUNCTION(JacobianDarcySystem2D)(void *ctx, CeedInt Q,
184184

185185
// -----------------------------------------------------------------------------
186186

187-
#endif //End of DARCY_MASS2D_H
187+
#endif //End of DARCY_SYSTEM2D_H

examples/Hdiv-mixed/qfunctions/darcy-system3d.h

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -195,4 +195,4 @@ CEED_QFUNCTION(JacobianDarcySystem3D)(void *ctx, CeedInt Q,
195195

196196
// -----------------------------------------------------------------------------
197197

198-
#endif //End of DARCY_MASS3D_H
198+
#endif //End of DARCY_SYSTEM3D_H
Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
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 RICHARD_FORCE2D_H
21+
#define RICHARD_FORCE2D_H
22+
23+
#include <math.h>
24+
#include "utils.h"
25+
26+
// See Matthew Farthing, Christopher Kees, Cass Miller (2003)
27+
// https://www.sciencedirect.com/science/article/pii/S0309170802001872
28+
// -----------------------------------------------------------------------------
29+
// Strong form:
30+
// k*K^{-1} * u = -\grad(p) + rho*g in \Omega x [0,T]
31+
// -\div(u) = -f + d (rho/rho_0*theta)/dt in \Omega x [0,T]
32+
// p = p_b on \Gamma_D x [0,T]
33+
// u.n = u_b on \Gamma_N x [0,T]
34+
// p = p_0 in \Omega, t = 0
35+
//
36+
// Where g is gravity vector, rho = rho_0*exp(beta * (p - p0)), p0 = 101325 Pa is atmospheric pressure
37+
// f = fs/rho_0, where g is gravity, rho_0 is the density at p_0, K = kappa*I, and
38+
// k_r = b_a + alpha_a * (\psi - x2), where \psi = p / (rho_0 * norm(g)) and x2 is vertical axis
39+
//
40+
// Weak form: Find (u, p) \in VxQ (V=H(div), Q=L^2) on \Omega
41+
// (v, k*K^{-1} * u) -(v, rho*g) - (\div(v), p) = - <v, p_b*n>_{\Gamma_D}
42+
// -(q, \div(u)) = -(q, f) + (v, d (rho/rho_0*theta)/dt )
43+
//
44+
// where k*K^{-1} = (rho_0^2*norm(g)/rho*k_r)*K^{-1}
45+
// This QFunction sets up the force
46+
// Inputs:
47+
// x : interpolation of the physical coordinate
48+
// w : weight of quadrature
49+
// J : dx/dX. x physical coordinate, X reference coordinate [-1,1]^dim
50+
//
51+
// Output:
52+
// force_u : which is 0.0 for this problem (-<v, p0 n> is in pressure-boundary qfunction)
53+
// force_p : -(q, f) = -\int( q * f * w*detJ)dx
54+
// -----------------------------------------------------------------------------
55+
// We have 3 experiment parameters as described in Table 1:P1, P2, P3
56+
// Matthew Farthing, Christopher Kees, Cass Miller (2003)
57+
// https://www.sciencedirect.com/science/article/pii/S0309170802001872
58+
#ifndef RICHARD_CTX
59+
#define RICHARD_CTX
60+
typedef struct RICHARDContext_ *RICHARDContext;
61+
struct RICHARDContext_ {
62+
CeedScalar kappa;
63+
CeedScalar alpha_a;
64+
CeedScalar b_a;
65+
CeedScalar rho_0;
66+
CeedScalar beta;
67+
CeedScalar g;
68+
CeedScalar p0;
69+
};
70+
#endif
71+
// -----------------------------------------------------------------------------
72+
// Force evaluation for Richard problem
73+
// -----------------------------------------------------------------------------
74+
CEED_QFUNCTION(RichardForce2D)(void *ctx, const CeedInt Q,
75+
const CeedScalar *const *in,
76+
CeedScalar *const *out) {
77+
// *INDENT-OFF*
78+
// Inputs
79+
const CeedScalar (*coords) = in[0],
80+
(*w) = in[1],
81+
(*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[2];
82+
// Outputs
83+
CeedScalar (*rhs_u) = out[0], (*rhs_p) = out[1],
84+
(*true_soln) = out[2];
85+
// Context
86+
RICHARDContext context = (RICHARDContext)ctx;
87+
const CeedScalar kappa = context->kappa;//10.;
88+
// Quadrature Point Loop
89+
CeedPragmaSIMD
90+
for (CeedInt i=0; i<Q; i++) {
91+
// Setup, (x,y) and J = dx/dX
92+
CeedScalar x = coords[i+0*Q], y = coords[i+1*Q];
93+
const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]},
94+
{dxdX[0][1][i], dxdX[1][1][i]}};
95+
const CeedScalar det_J = MatDet2x2(J);
96+
// *INDENT-ON*
97+
CeedScalar pe = sin(PI_DOUBLE*x) * sin(PI_DOUBLE*y);
98+
CeedScalar grad_pe[2] = {PI_DOUBLE*cos(PI_DOUBLE*x) *sin(PI_DOUBLE*y), PI_DOUBLE*sin(PI_DOUBLE*x) *cos(PI_DOUBLE*y)};
99+
CeedScalar K[2][2] = {{kappa, 0.},{0., kappa}};
100+
CeedScalar ue[2];
101+
AlphaMatVecMult2x2(-1., K, grad_pe, ue);
102+
CeedScalar f = 2*PI_DOUBLE*PI_DOUBLE*sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
103+
104+
// 1st eq: component 1
105+
rhs_u[i+0*Q] = 0.;
106+
// 1st eq: component 2
107+
rhs_u[i+1*Q] = 0.;
108+
// 2nd eq
109+
rhs_p[i] = -f*w[i]*det_J;
110+
// True solution Ue=[p,u]
111+
true_soln[i+0*Q] = pe;
112+
true_soln[i+1*Q] = ue[0];
113+
true_soln[i+2*Q] = ue[1];
114+
} // End of Quadrature Point Loop
115+
return 0;
116+
}
117+
// -----------------------------------------------------------------------------
118+
119+
#endif //End of RICHARD_FORCE2D_H

0 commit comments

Comments
 (0)