Skip to content

Commit 833e410

Browse files
rezgarshakerirezgarshakeri
committed
Setup Richard true solution
1 parent 5c2906a commit 833e410

7 files changed

Lines changed: 92 additions & 97 deletions

File tree

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ struct AppCtx_ {
1818
// Problem type arguments
1919
PetscFunctionList problems;
2020
char problem_name[PETSC_MAX_PATH_LEN];
21+
CeedContextFieldLabel solution_time_label;
2122
};
2223

2324
// 2) richard

‎examples/Hdiv-mixed/main.c‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -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
@@ -143,6 +143,7 @@ int main(int argc, char **argv) {
143143
// ---------------------------------------------------------------------------
144144
Vec U0;
145145
CreateInitialConditions(dm, ceed_data, &U0);
146+
VecView(U0, PETSC_VIEWER_STDOUT_WORLD);
146147
PetscCall( VecDestroy(&U0) );
147148
}
148149

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

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@
1919

2020
#include "../include/register-problem.h"
2121
#include "../qfunctions/richard-system2d.h"
22-
#include "../qfunctions/richard-mms2d.h"
22+
#include "../qfunctions/richard-true2d.h"
2323
#include "../qfunctions/pressure-boundary2d.h"
2424
#include "petscsystypes.h"
2525

@@ -41,8 +41,8 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
4141
problem_data->quadrature_mode = CEED_GAUSS;
4242
problem_data->ics = RichardICs2D;
4343
problem_data->ics_loc = RichardICs2D_loc;
44-
//problem_data->force = RichardTrue2D;
45-
//problem_data->force_loc = RichardTrue2D_loc;
44+
problem_data->true_solution = RichardTrue2D;
45+
problem_data->true_solution_loc = RichardTrue2D_loc;
4646
//problem_data->residual = RichardSystem2D;
4747
//problem_data->residual_loc = RichardSystem2D_loc;
4848
//problem_data->jacobian = JacobianRichardSystem2D;
@@ -81,11 +81,12 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) {
8181
richard_ctx->beta = beta;
8282
richard_ctx->g = g;
8383
richard_ctx->p0 = p0;
84-
richard_ctx->time = 0.;
8584
richard_ctx->gamma = 5.;
8685
CeedQFunctionContextCreate(ceed, &richard_context);
8786
CeedQFunctionContextSetData(richard_context, CEED_MEM_HOST, CEED_COPY_VALUES,
8887
sizeof(*richard_ctx), richard_ctx);
88+
CeedQFunctionContextRegisterDouble(richard_context, "time",
89+
offsetof(struct RICHARDContext_, t), 1, "current solver time");
8990
problem_data->qfunction_context = richard_context;
9091
CeedQFunctionContextSetDataDestroy(richard_context, CEED_MEM_HOST,
9192
FreeContextPetsc);

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

Lines changed: 4 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -67,13 +67,11 @@
6767
typedef struct RICHARDContext_ *RICHARDContext;
6868
struct RICHARDContext_ {
6969
CeedScalar kappa;
70-
CeedScalar alpha_a;
71-
CeedScalar b_a;
72-
CeedScalar rho_a0;
73-
CeedScalar beta;
7470
CeedScalar g;
75-
CeedScalar p0;
76-
CeedScalar time;
71+
CeedScalar rho_a0;
72+
CeedScalar alpha_a, b_a;
73+
CeedScalar beta, p0;
74+
CeedScalar t;
7775
CeedScalar gamma;
7876
};
7977
#endif

examples/Hdiv-mixed/qfunctions/richard-mms2d.h renamed to examples/Hdiv-mixed/qfunctions/richard-true2d.h

Lines changed: 34 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -17,19 +17,18 @@
1717
/// @file
1818
/// Force of Richard problem 2D (quad element) using PETSc
1919

20-
#ifndef RICHARD_MMS2D_H
21-
#define RICHARD_MMS2D_H
20+
#ifndef RICHARD_TRUE2D_H
21+
#define RICHARD_TRUE2D_H
2222

2323
#include <math.h>
24-
#include "ceed/ceed-f64.h"
2524
#include "utils.h"
2625

2726
// See Matthew Farthing, Christopher Kees, Cass Miller (2003)
2827
// https://www.sciencedirect.com/science/article/pii/S0309170802001872
2928
// -----------------------------------------------------------------------------
3029
// Strong form:
3130
// u = -rho*k_r*K *[grad(\psi) - rho*g_u] in \Omega x [0,T]
32-
// -\div(u) = -f + d (rho*\theta)/dt in \Omega x [0,T]
31+
// -\div(u) = -f + d (rho*theta)/dt in \Omega x [0,T]
3332
// p = p_b on \Gamma_D x [0,T]
3433
// u.n = u_b on \Gamma_N x [0,T]
3534
// p = p_0 in \Omega, t = 0
@@ -45,15 +44,13 @@
4544
// We solve MMS for K = kappa*I and beta=0 ==> rho=1 and \theta = alpha_a*\psi, so
4645
// -(q, d (rho*\theta)/dt ) = -alpha_a*(q, d(\psi)/dt )
4746
//
48-
// This QFunction sets up the force
47+
// This QFunction setup the true solution and forcing f of the above equation
4948
// Inputs:
50-
// x : interpolation of the physical coordinate
51-
// w : weight of quadrature
52-
// J : dx/dX. x physical coordinate, X reference coordinate [-1,1]^dim
49+
// coords: physical coordinate
5350
//
5451
// Output:
55-
// force_u : which is 0.0 for this problem (-<v, p0 n> is in pressure-boundary qfunction)
56-
// force_p : 0.0
52+
// true_force : = div(u) + d (rho*theta)/dt
53+
// true_solution : = [\psi, u] where \psi, u are the exact solution solution
5754
// -----------------------------------------------------------------------------
5855
// We have 3 experiment parameters as described in Table 1:P1, P2, P3
5956
// Matthew Farthing, Christopher Kees, Cass Miller (2003)
@@ -63,13 +60,11 @@
6360
typedef struct RICHARDContext_ *RICHARDContext;
6461
struct RICHARDContext_ {
6562
CeedScalar kappa;
66-
CeedScalar alpha_a;
67-
CeedScalar b_a;
68-
CeedScalar rho_a0;
69-
CeedScalar beta;
7063
CeedScalar g;
71-
CeedScalar p0;
72-
CeedScalar time;
64+
CeedScalar rho_a0;
65+
CeedScalar alpha_a, b_a;
66+
CeedScalar beta, p0;
67+
CeedScalar t;
7368
CeedScalar gamma;
7469
};
7570
#endif
@@ -99,7 +94,7 @@ CEED_QFUNCTION(RichardICs2D)(void *ctx, const CeedInt Q,
9994
} // End of Quadrature Point Loop
10095
return 0;
10196
}
102-
/*
97+
10398
// -----------------------------------------------------------------------------
10499
// True solution for Richard problem
105100
// -----------------------------------------------------------------------------
@@ -108,50 +103,43 @@ CEED_QFUNCTION(RichardTrue2D)(void *ctx, const CeedInt Q,
108103
CeedScalar *const *out) {
109104
// *INDENT-OFF*
110105
// Inputs
111-
const CeedScalar (*coords) = in[0],
112-
(*w) = in[1],
113-
(*dxdX)[2][CEED_Q_VLA] = (const CeedScalar(*)[2][CEED_Q_VLA])in[2];
106+
const CeedScalar (*coords) = in[0];
114107
// Outputs
115-
CeedScalar (*true_soln) = out[0];
108+
CeedScalar (*true_force) = out[0], (*true_solution) = out[1];
116109
// Context
117110
RICHARDContext context = (RICHARDContext)ctx;
111+
const CeedScalar kappa = context->kappa;
112+
const CeedScalar alpha_a = context->alpha_a;
113+
const CeedScalar b_a = context->b_a;
118114
const CeedScalar gamma = context->gamma;
119-
const CeedScalar rho_a0 = context->rho_a0;
120-
const CeedScalar g = context->g;
121-
CeedScalar t = context->time;
115+
CeedScalar t = context->t;
116+
printf("time %f \n", t);
122117
// Quadrature Point Loop
123118
CeedPragmaSIMD
124119
for (CeedInt i=0; i<Q; i++) {
125-
// Setup, (x,y) and J = dx/dX
126120
CeedScalar x = coords[i+0*Q], y = coords[i+1*Q];
127-
const CeedScalar J[2][2] = {{dxdX[0][0][i], dxdX[1][0][i]},
128-
{dxdX[0][1][i], dxdX[1][1][i]}};
129-
const CeedScalar det_J = MatDet2x2(J);
130-
131121
CeedScalar psi = exp(-gamma*t)*sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
132122
CeedScalar psi_x = PI_DOUBLE*exp(-gamma*t)*cos(PI_DOUBLE*x)*sin(PI_DOUBLE*y);
133-
CeedScalar psi_xx = -PI_DOUBLE*PI_DOUBLE*psi;
134123
CeedScalar psi_y = PI_DOUBLE*exp(-gamma*t)*sin(PI_DOUBLE*x)*cos(PI_DOUBLE*y);
135-
CeedScalar psi_yy = -PI_DOUBLE*PI_DOUBLE*psi;
136-
137-
// k_r = b_a + alpha_a * (psi - x2)
138-
CeedScalar k_r = b_a + alpha_a * (psi - y);
139-
CeedScalar u[2] = {-k_r*kappa*psi_x, -k_r*kappa*(psi_y-1)};
140-
CeedScalar div_u = -kappa*(alpha_a*psi_x*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi
141-
+ alpha_a*(psi_y-1)*(psi_y-1)
142-
- k_r*PI_DOUBLE*PI_DOUBLE*psi);
143124

144-
CeedScalar f = div_u -alpha_a*gamma*psi;
145-
AlphaMatVecMult2x2(-1./k, K, grad_pe, ue);
125+
// k_r = b_a + alpha_a * (1 - x*y)
126+
CeedScalar k_r = b_a + alpha_a*(1-x*y);
127+
// rho = rho_a/rho_a0
128+
CeedScalar rho = 1.;
129+
// u = -rho*k_r*K *[grad(\psi) - rho*g_u]
130+
CeedScalar u[2] = {-rho*k_r*kappa*psi_x, -rho*k_r*kappa*(psi_y-1)};
131+
CeedScalar div_u = -rho*kappa*(-alpha_a*y*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi
132+
-alpha_a*x*psi_y - k_r*PI_DOUBLE*PI_DOUBLE*psi);
146133

147-
// True solution Ue=[p,u]
148-
true_soln[i+0*Q] = psi;
149-
true_soln[i+1*Q] = u[0];
150-
true_soln[i+2*Q] = u[1];
134+
// True Force: f = \div(u) + d (rho*theta)/dt
135+
true_force[i+0*Q] = div_u -alpha_a*gamma*psi;
136+
// True Solution
137+
true_solution[i+0*Q] = psi;
138+
true_solution[i+1*Q] = u[0];
139+
true_solution[i+2*Q] = u[1];
151140
} // End of Quadrature Point Loop
152141
return 0;
153142
}
154-
*/
155143
// -----------------------------------------------------------------------------
156144

157-
#endif //End of RICHARD_MMS2D_H
145+
#endif //End of RICHARD_TRUE2D_H

‎examples/Hdiv-mixed/src/setup-libceed.c‎

Lines changed: 35 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -245,33 +245,6 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
245245
(PetscScalar *)coordArray);
246246
PetscCall( VecRestoreArrayRead(coords, &coordArray) );
247247

248-
// ---------------------------------------------------------------------------
249-
// Setup true solution and force
250-
// ---------------------------------------------------------------------------
251-
CeedVector true_vec, true_force;
252-
CeedVectorCreate(ceed, num_elem*num_qpts*(dim+1), &true_vec);
253-
CeedVectorCreate(ceed, num_elem*num_qpts*1, &true_force);
254-
// Create the q-function that sets up the RHS and true solution
255-
CeedQFunctionCreateInterior(ceed, 1, problem_data->true_solution,
256-
problem_data->true_solution_loc, &qf_true);
257-
CeedQFunctionSetContext(qf_true, problem_data->qfunction_context);
258-
//CeedQFunctionContextDestroy(&problem_data->qfunction_context);
259-
CeedQFunctionAddInput(qf_true, "x", num_comp_x, CEED_EVAL_INTERP);
260-
CeedQFunctionAddOutput(qf_true, "true force", 1, CEED_EVAL_NONE);
261-
CeedQFunctionAddOutput(qf_true, "true solution", dim+1, CEED_EVAL_NONE);
262-
// Create the operator that builds the RHS and true solution
263-
CeedOperatorCreate(ceed, qf_true, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
264-
&op_true);
265-
CeedOperatorSetField(op_true, "x", ceed_data->elem_restr_x,
266-
ceed_data->basis_x, ceed_data->x_coord);
267-
CeedOperatorSetField(op_true, "true force", ceed_data->elem_restr_p_i,
268-
CEED_BASIS_COLLOCATED, true_force);
269-
CeedOperatorSetField(op_true, "true solution", ceed_data->elem_restr_U_i,
270-
CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
271-
// Setup RHS and true solution
272-
CeedOperatorApply(op_true, ceed_data->x_coord, true_vec,
273-
CEED_REQUEST_IMMEDIATE);
274-
275248
if (problem_data->has_ts) {
276249
// ---------------------------------------------------------------------------
277250
// Setup qfunction for initial conditions
@@ -299,12 +272,46 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
299272
ceed_data->qf_ics = qf_ics;
300273
ceed_data->op_ics = op_ics;
301274
}
275+
276+
// ---------------------------------------------------------------------------
277+
// Setup true solution and force
278+
// ---------------------------------------------------------------------------
279+
CeedVector true_vec, true_force;
280+
CeedVectorCreate(ceed, num_elem*num_qpts*(dim+1), &true_vec);
281+
CeedVectorCreate(ceed, num_elem*num_qpts*1, &true_force);
282+
// Create the q-function that sets up the RHS and true solution
283+
CeedQFunctionCreateInterior(ceed, 1, problem_data->true_solution,
284+
problem_data->true_solution_loc, &qf_true);
285+
CeedQFunctionSetContext(qf_true, problem_data->qfunction_context);
286+
CeedQFunctionAddInput(qf_true, "x", num_comp_x, CEED_EVAL_INTERP);
287+
CeedQFunctionAddOutput(qf_true, "true force", 1, CEED_EVAL_NONE);
288+
CeedQFunctionAddOutput(qf_true, "true solution", dim+1, CEED_EVAL_NONE);
289+
// Create the operator that builds the RHS and true solution
290+
CeedOperatorCreate(ceed, qf_true, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
291+
&op_true);
292+
if (problem_data->has_ts) {
293+
double value_time = 2.0;
294+
CeedOperatorContextGetFieldLabel(op_true, "time",
295+
&app_ctx->solution_time_label);
296+
CeedOperatorContextSetDouble(op_true,
297+
app_ctx->solution_time_label, &value_time);
298+
}
299+
CeedOperatorSetField(op_true, "x", ceed_data->elem_restr_x,
300+
ceed_data->basis_x, ceed_data->x_coord);
301+
CeedOperatorSetField(op_true, "true force", ceed_data->elem_restr_p_i,
302+
CEED_BASIS_COLLOCATED, true_force);
303+
CeedOperatorSetField(op_true, "true solution", ceed_data->elem_restr_U_i,
304+
CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
305+
// Setup RHS and true solution
306+
CeedOperatorApply(op_true, ceed_data->x_coord, true_vec,
307+
CEED_REQUEST_IMMEDIATE);
302308
// ---------------------------------------------------------------------------
303309
// Persistent libCEED vectors
304310
// ---------------------------------------------------------------------------
305311
// -- Operator action variables: we use them in setup-solvers.c
306312
CeedVectorCreate(ceed, U_loc_size, &ceed_data->x_ceed);
307313
CeedVectorCreate(ceed, U_loc_size, &ceed_data->y_ceed);
314+
if (!problem_data->has_ts) {
308315
// Local residual evaluator
309316
// ---------------------------------------------------------------------------
310317
// Create the QFunction and Operator that computes the residual of the PDE.
@@ -443,6 +450,7 @@ PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
443450
CeedVectorDestroy(&true_force);
444451
CeedQFunctionDestroy(&qf_true);
445452
CeedOperatorDestroy(&op_true);
453+
}
446454
PetscFunctionReturn(0);
447455
};
448456
// -----------------------------------------------------------------------------

‎examples/Hdiv-mixed/src/setup-ts.c‎

Lines changed: 11 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *U0) {
4444

4545
}
4646
// -----------------------------------------------------------------------------
47-
// Setup operator context data
47+
// Setup operator context data for Richard problem
4848
// -----------------------------------------------------------------------------
4949
PetscErrorCode SetupResidualOperatorCtx_Ut(DM dm, Ceed ceed, CeedData ceed_data,
5050
OperatorApplyContext ctx_residual_ut) {
@@ -68,22 +68,20 @@ PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
6868
OperatorApplyContext ctx = (OperatorApplyContext)ctx_residual_ut;
6969
const PetscScalar *x, *x_t;
7070
PetscScalar *y;
71-
Vec X_loc = ctx->X_loc, X_t_loc = ctx->X_t_loc,
72-
Y_loc = ctx->Y_loc;
7371
PetscMemType x_mem_type, x_t_mem_type, y_mem_type;
7472
PetscFunctionBeginUser;
7573

7674
// Update time dependent data
7775

7876

7977
// Global-to-local
80-
PetscCall( DMGlobalToLocal(ctx->dm, X, INSERT_VALUES, X_loc) );
81-
PetscCall( DMGlobalToLocal(ctx->dm, X_t, INSERT_VALUES, X_t_loc) );
78+
PetscCall( DMGlobalToLocal(ctx->dm, X, INSERT_VALUES, ctx->X_loc) );
79+
PetscCall( DMGlobalToLocal(ctx->dm, X_t, INSERT_VALUES, ctx->X_t_loc) );
8280

8381
// Place PETSc vectors in CEED vectors
84-
PetscCall( VecGetArrayReadAndMemType(X_loc, &x, &x_mem_type) );
85-
PetscCall( VecGetArrayReadAndMemType(X_t_loc, &x_t, &x_t_mem_type) );
86-
PetscCall( VecGetArrayAndMemType(Y_loc, &y, &y_mem_type) );
82+
PetscCall( VecGetArrayReadAndMemType(ctx->X_loc, &x, &x_mem_type) );
83+
PetscCall( VecGetArrayReadAndMemType(ctx->X_t_loc, &x_t, &x_t_mem_type) );
84+
PetscCall( VecGetArrayAndMemType(ctx->Y_loc, &y, &y_mem_type) );
8785
CeedVectorSetArray(ctx->x_ceed, MemTypeP2C(x_mem_type), CEED_USE_POINTER,
8886
(PetscScalar *)x);
8987
CeedVectorSetArray(ctx->x_t_ceed, MemTypeP2C(x_t_mem_type),
@@ -98,16 +96,16 @@ PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
9896
CeedVectorTakeArray(ctx->x_ceed, MemTypeP2C(x_mem_type), NULL);
9997
CeedVectorTakeArray(ctx->x_t_ceed, MemTypeP2C(x_t_mem_type), NULL);
10098
CeedVectorTakeArray(ctx->y_ceed, MemTypeP2C(y_mem_type), NULL);
101-
PetscCall( VecRestoreArrayReadAndMemType(X_loc, &x) );
102-
PetscCall( VecRestoreArrayReadAndMemType(X_t_loc, &x_t) );
103-
PetscCall( VecRestoreArrayAndMemType(Y_loc, &y) );
99+
PetscCall( VecRestoreArrayReadAndMemType(ctx->X_loc, &x) );
100+
PetscCall( VecRestoreArrayReadAndMemType(ctx->X_t_loc, &x_t) );
101+
PetscCall( VecRestoreArrayAndMemType(ctx->Y_loc, &y) );
104102

105103
// Local-to-Global
106104
PetscCall( VecZeroEntries(Y) );
107-
PetscCall( DMLocalToGlobal(ctx->dm, Y_loc, ADD_VALUES, Y) );
105+
PetscCall( DMLocalToGlobal(ctx->dm, ctx->Y_loc, ADD_VALUES, Y) );
108106

109107
// Restore vectors
110-
PetscCall( DMRestoreLocalVector(ctx->dm, &Y_loc) );
108+
PetscCall( DMRestoreLocalVector(ctx->dm, &ctx->Y_loc) );
111109

112110
PetscFunctionReturn(0);
113111
}

0 commit comments

Comments
 (0)