Skip to content

Commit 688e27d

Browse files
rezgarshakerirezgarshakeri
committed
Moved everything to left and deleted qfunction for rhs
1 parent 5a5d3c0 commit 688e27d

28 files changed

+1156
-567
lines changed

‎examples/Hdiv-mixed/conv_test_result.csv

Lines changed: 11 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
11
run,mesh_res,error_u,error_p
2-
0,2,0.85636,0.09544
3-
1,3,0.43484,0.05265
4-
2,4,0.26359,0.03280
5-
3,5,0.17618,0.02214
6-
4,6,0.12591,0.01588
7-
5,7,0.09423,0.01191
8-
6,8,0.07307,0.00925
9-
7,9,0.05826,0.00738
10-
8,10,0.04750,0.00602
11-
9,11,0.03946,0.00501
12-
10,12,0.03328,0.00422
2+
0,2,9.13818,0.09469
3+
1,3,4.60103,0.05186
4+
2,4,2.75883,0.03199
5+
3,5,1.81813,0.02132
6+
4,6,1.27780,0.01505
7+
5,7,0.93759,0.01108
8+
6,8,0.71073,0.00842
9+
7,9,0.55232,0.00655
10+
8,10,0.43767,0.00520
11+
9,11,0.35229,0.00418
12+
10,12,0.28723,0.00340

‎examples/Hdiv-mixed/convrate_mixed.png

-1.57 KB
Loading

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

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx);
1717
// 3) darcy3dprism
1818

1919
// 4) richard
20+
PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx);
2021

2122
extern int FreeContextPetsc(void *);
2223

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

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,5 @@ PetscErrorCode CreateRestrictionFromPlexOriented(Ceed ceed, DM dm, CeedInt P,
1919
// Set up libCEED for a given degree
2020
PetscErrorCode SetupLibceed(DM dm, Ceed ceed, AppCtx app_ctx,
2121
ProblemData problem_data,
22-
PetscInt rhs_loc_size, CeedData ceed_data,
23-
CeedVector rhs_ceed, CeedVector *target);
22+
PetscInt U_loc_size, CeedData ceed_data);
2423
#endif // setuplibceed_h

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

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -16,10 +16,8 @@ PetscErrorCode ApplyJacobian(Mat A, Vec X, Vec Y);
1616
PetscErrorCode SNESFormResidual(SNES snes, Vec X, Vec Y, void *ctx);
1717
PetscErrorCode SNESFormJacobian(SNES snes, Vec U, Mat J, Mat J_pre, void *ctx);
1818
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);
19+
VecType vec_type, SNES snes, KSP ksp, Vec *U);
2120
PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U,
22-
CeedVector target,
2321
CeedScalar *l2_error_u, CeedScalar *l2_error_p);
2422
PetscErrorCode PrintOutput(Ceed ceed,
2523
CeedMemType mem_type_backend,

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

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
#ifndef setup_ts_h
2+
#define setup_ts_h
3+
4+
#include <ceed.h>
5+
#include <petsc.h>
6+
7+
#include "structs.h"
8+
9+
PetscErrorCode CreateInitialConditions(DM dm, CeedData ceed_data, Vec *U0);
10+
PetscErrorCode SetupResidualOperatorCtx_Ut(DM dm, Ceed ceed, CeedData ceed_data,
11+
OperatorApplyContext ctx_residual_ut);
12+
PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y,
13+
void *ctx_residual_ut);
14+
PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx,
15+
Vec *U, PetscScalar *f_time, TS *ts);
16+
17+
#endif // setup_ts_h

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

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -67,8 +67,8 @@ struct Physics_ {
6767
typedef struct OperatorApplyContext_ *OperatorApplyContext;
6868
struct OperatorApplyContext_ {
6969
MPI_Comm comm;
70-
Vec X_loc, Y_loc;
71-
CeedVector x_ceed, y_ceed;
70+
Vec X_loc, Y_loc, X_t_loc;
71+
CeedVector x_ceed, y_ceed, x_t_ceed;
7272
CeedOperator op_apply;
7373
DM dm;
7474
Ceed ceed;
@@ -79,23 +79,24 @@ typedef struct CeedData_ *CeedData;
7979
struct CeedData_ {
8080
CeedBasis basis_x, basis_u, basis_p, basis_u_face;
8181
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;
82+
elem_restr_p, elem_restr_p_i;
83+
CeedQFunction qf_residual, qf_jacobian, qf_error, qf_ics;
84+
CeedOperator op_residual, op_jacobian, op_error, op_ics;
85+
CeedVector x_ceed, y_ceed, x_coord, U0_ceed, x_t_ceed;
86+
OperatorApplyContext ctx_residual, ctx_jacobian, ctx_error, ctx_residual_ut;
8787
};
8888

8989
// Problem specific data
9090
typedef struct ProblemData_ *ProblemData;
9191
struct ProblemData_ {
92-
CeedQFunctionUser force, residual, jacobian, error,
93-
setup_true, bc_pressure;
94-
const char *force_loc, *residual_loc, *jacobian_loc,
95-
*error_loc, *setup_true_loc, *bc_pressure_loc;
92+
CeedQFunctionUser true_solution, residual, jacobian, error, ics,
93+
bc_pressure;
94+
const char *true_solution_loc, *residual_loc, *jacobian_loc,
95+
*error_loc, *bc_pressure_loc, *ics_loc;
9696
CeedQuadMode quadrature_mode;
9797
CeedInt elem_node, dim, q_data_size_face;
9898
CeedQFunctionContext qfunction_context;
99+
PetscBool has_ts;
99100
};
100101

101102
#endif // structs_h

‎examples/Hdiv-mixed/main.c

Lines changed: 25 additions & 38 deletions
Original file line numberDiff line numberDiff line change
@@ -76,9 +76,6 @@ int main(int argc, char **argv) {
7676
Physics phys_ctx;
7777
PetscCall( PetscCalloc1(1, &phys_ctx) );
7878

79-
OperatorApplyContext op_apply_ctx;
80-
PetscCall( PetscCalloc1(1, &op_apply_ctx) );
81-
8279
// ---------------------------------------------------------------------------
8380
// Process command line options
8481
// ---------------------------------------------------------------------------
@@ -119,39 +116,35 @@ int main(int argc, char **argv) {
119116
// ---------------------------------------------------------------------------
120117
// Create local Force vector
121118
// ---------------------------------------------------------------------------
122-
Vec F_loc;
123-
PetscInt F_loc_size;
124-
CeedScalar *f;
125-
CeedVector force_ceed, target, bc_pressure;
126-
PetscMemType force_mem_type;
127-
PetscCall( DMCreateLocalVector(dm, &F_loc) );
119+
Vec U_loc;
120+
PetscInt U_loc_size;
121+
//CeedVector bc_pressure;
122+
PetscCall( DMCreateLocalVector(dm, &U_loc) );
128123
// Local size for libCEED
129-
PetscCall( VecGetSize(F_loc, &F_loc_size) );
130-
PetscCall( VecZeroEntries(F_loc) );
131-
PetscCall( VecGetArrayAndMemType(F_loc, &f, &force_mem_type) );
132-
CeedVectorCreate(ceed, F_loc_size, &force_ceed);
133-
CeedVectorSetArray(force_ceed, MemTypeP2C(force_mem_type), CEED_USE_POINTER, f);
134-
CeedVectorCreate(ceed, F_loc_size, &bc_pressure);
135-
CeedVectorSetArray(bc_pressure, MemTypeP2C(force_mem_type), CEED_USE_POINTER,
136-
f);
124+
PetscCall( VecGetSize(U_loc, &U_loc_size) );
125+
137126
// ---------------------------------------------------------------------------
138-
// Setup libCEED - Compute local F and true solution (target)
127+
// Setup libCEED
139128
// ---------------------------------------------------------------------------
140129
// -- Set up libCEED objects
141130
PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
142-
F_loc_size, ceed_data, force_ceed, &target) );
131+
U_loc_size, ceed_data) );
143132
//CeedVectorView(force_ceed, "%12.8f", stdout);
144-
PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
145-
bc_pressure) );
133+
//PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
134+
// bc_pressure) );
135+
136+
146137
// ---------------------------------------------------------------------------
147-
// Create global F
138+
// Setup TSSolve for Richard problem
148139
// ---------------------------------------------------------------------------
149-
Vec F;
150-
CeedVectorTakeArray(force_ceed, MemTypeP2C(force_mem_type), NULL);
151-
PetscCall( VecRestoreArrayAndMemType(F_loc, &f) );
152-
PetscCall( DMCreateGlobalVector(dm, &F) );
153-
PetscCall( VecZeroEntries(F) );
154-
PetscCall( DMLocalToGlobal(dm, F_loc, ADD_VALUES, F) );
140+
if (problem_data->has_ts) {
141+
// ---------------------------------------------------------------------------
142+
// Create global initial conditions
143+
// ---------------------------------------------------------------------------
144+
Vec U0;
145+
CreateInitialConditions(dm, ceed_data, &U0);
146+
PetscCall( VecDestroy(&U0) );
147+
}
155148

156149
// ---------------------------------------------------------------------------
157150
// Solve PDE
@@ -162,14 +155,14 @@ int main(int argc, char **argv) {
162155
Vec U;
163156
PetscCall( SNESCreate(comm, &snes) );
164157
PetscCall( SNESGetKSP(snes, &ksp) );
165-
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, F, &U) );
158+
PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) );
166159
//VecView(U, PETSC_VIEWER_STDOUT_WORLD);
167160

168161
// ---------------------------------------------------------------------------
169162
// Compute L2 error of mms problem
170163
// ---------------------------------------------------------------------------
171164
CeedScalar l2_error_u, l2_error_p;
172-
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, target, &l2_error_u,
165+
PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u,
173166
&l2_error_p) );
174167

175168
// ---------------------------------------------------------------------------
@@ -194,11 +187,8 @@ int main(int argc, char **argv) {
194187
// Free PETSc objects
195188
PetscCall( DMDestroy(&dm) );
196189
PetscCall( VecDestroy(&U) );
197-
PetscCall( VecDestroy(&F) );
198-
PetscCall( VecDestroy(&F_loc) );
190+
PetscCall( VecDestroy(&U_loc) );
199191
PetscCall( SNESDestroy(&snes) );
200-
PetscCall( VecDestroy(&op_apply_ctx->Y_loc) );
201-
PetscCall( VecDestroy(&op_apply_ctx->X_loc) );
202192

203193
// -- Function list
204194
PetscCall( PetscFunctionListDestroy(&app_ctx->problems) );
@@ -207,12 +197,9 @@ int main(int argc, char **argv) {
207197
PetscCall( PetscFree(app_ctx) );
208198
PetscCall( PetscFree(problem_data) );
209199
PetscCall( PetscFree(phys_ctx) );
210-
PetscCall( PetscFree(op_apply_ctx) );
211200

212201
// Free libCEED objects
213-
CeedVectorDestroy(&force_ceed);
214-
CeedVectorDestroy(&bc_pressure);
215-
CeedVectorDestroy(&target);
202+
//CeedVectorDestroy(&bc_pressure);
216203
PetscCall( CeedDataDestroy(ceed_data) );
217204
CeedDestroy(&ceed);
218205

‎examples/Hdiv-mixed/main.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,5 +10,6 @@
1010
#include "include/register-problem.h"
1111
#include "include/setup-matops.h"
1212
#include "include/setup-solvers.h"
13+
#include "include/setup-ts.h"
1314

1415
#endif // MAIN_H

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

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
/// Utility functions for setting up Darcy problem in 2D
1919

2020
#include "../include/register-problem.h"
21-
#include "../qfunctions/darcy-force2d.h"
21+
#include "../qfunctions/darcy-true2d.h"
2222
#include "../qfunctions/darcy-system2d.h"
2323
#include "../qfunctions/darcy-error2d.h"
2424
#include "../qfunctions/pressure-boundary2d.h"
@@ -39,8 +39,8 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
3939
problem_data->elem_node = 4;
4040
problem_data->q_data_size_face = 3;
4141
problem_data->quadrature_mode = CEED_GAUSS;
42-
problem_data->force = DarcyForce2D;
43-
problem_data->force_loc = DarcyForce2D_loc;
42+
problem_data->true_solution = DarcyTrue2D;
43+
problem_data->true_solution_loc = DarcyTrue2D_loc;
4444
problem_data->residual = DarcySystem2D;
4545
problem_data->residual_loc = DarcySystem2D_loc;
4646
problem_data->jacobian = JacobianDarcySystem2D;
@@ -49,23 +49,37 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) {
4949
problem_data->error_loc = DarcyError2D_loc;
5050
problem_data->bc_pressure = BCPressure2D;
5151
problem_data->bc_pressure_loc = BCPressure2D_loc;
52+
problem_data->has_ts = PETSC_FALSE;
5253

5354
// ------------------------------------------------------
5455
// Command line Options
5556
// ------------------------------------------------------
56-
CeedScalar kappa = 1.;
57+
CeedScalar kappa = 1., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.;
5758
PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
5859
PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
5960
kappa, &kappa, NULL));
61+
PetscCall( PetscOptionsScalar("-rho_a0", "Density at p0", NULL,
62+
rho_a0, &rho_a0, NULL));
63+
PetscCall( PetscOptionsScalar("-alpha_a", "Parameter for relative permeability",
64+
NULL,
65+
alpha_a, &alpha_a, NULL));
66+
PetscCall( PetscOptionsScalar("-b_a", "Parameter for relative permeability",
67+
NULL,
68+
b_a, &b_a, NULL));
6069
PetscOptionsEnd();
6170

6271
darcy_ctx->kappa = kappa;
72+
darcy_ctx->rho_a0 = rho_a0;
73+
darcy_ctx->g = g;
74+
darcy_ctx->alpha_a = alpha_a;
75+
darcy_ctx->b_a = b_a;
76+
6377
CeedQFunctionContextCreate(ceed, &darcy_context);
6478
CeedQFunctionContextSetData(darcy_context, CEED_MEM_HOST, CEED_COPY_VALUES,
6579
sizeof(*darcy_ctx), darcy_ctx);
6680
problem_data->qfunction_context = darcy_context;
6781
CeedQFunctionContextSetDataDestroy(darcy_context, CEED_MEM_HOST,
6882
FreeContextPetsc);
69-
//PetscCall( PetscFree(darcy_ctx) );
83+
PetscCall( PetscFree(darcy_ctx) );
7084
PetscFunctionReturn(0);
7185
}

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

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
/// Utility functions for setting up Darcy problem in 3D
1919

2020
#include "../include/register-problem.h"
21-
#include "../qfunctions/darcy-force3d.h"
21+
#include "../qfunctions/darcy-true3d.h"
2222
#include "../qfunctions/darcy-system3d.h"
2323
#include "../qfunctions/darcy-error3d.h"
2424
#include "../qfunctions/pressure-boundary3d.h"
@@ -39,8 +39,8 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx) {
3939
problem_data->elem_node = 8;
4040
problem_data->q_data_size_face = 4;
4141
problem_data->quadrature_mode = CEED_GAUSS;
42-
problem_data->force = DarcyForce3D;
43-
problem_data->force_loc = DarcyForce3D_loc;
42+
problem_data->true_solution = DarcyTrue3D;
43+
problem_data->true_solution_loc = DarcyTrue3D_loc;
4444
problem_data->residual = DarcySystem3D;
4545
problem_data->residual_loc = DarcySystem3D_loc;
4646
problem_data->jacobian = JacobianDarcySystem3D;
@@ -49,23 +49,37 @@ PetscErrorCode Hdiv_DARCY3D(Ceed ceed, ProblemData problem_data, void *ctx) {
4949
problem_data->error_loc = DarcyError3D_loc;
5050
problem_data->bc_pressure = BCPressure3D;
5151
problem_data->bc_pressure_loc = BCPressure3D_loc;
52+
problem_data->has_ts = PETSC_FALSE;
5253

5354
// ------------------------------------------------------
5455
// Command line Options
5556
// ------------------------------------------------------
56-
CeedScalar kappa = 1.;
57+
CeedScalar kappa = 1., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.;
5758
PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL);
5859
PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL,
5960
kappa, &kappa, NULL));
61+
PetscCall( PetscOptionsScalar("-rho_a0", "Density at p0", NULL,
62+
rho_a0, &rho_a0, NULL));
63+
PetscCall( PetscOptionsScalar("-alpha_a", "Parameter for relative permeability",
64+
NULL,
65+
alpha_a, &alpha_a, NULL));
66+
PetscCall( PetscOptionsScalar("-b_a", "Parameter for relative permeability",
67+
NULL,
68+
b_a, &b_a, NULL));
6069
PetscOptionsEnd();
6170

6271
darcy_ctx->kappa = kappa;
72+
darcy_ctx->rho_a0 = rho_a0;
73+
darcy_ctx->g = g;
74+
darcy_ctx->alpha_a = alpha_a;
75+
darcy_ctx->b_a = b_a;
76+
6377
CeedQFunctionContextCreate(ceed, &darcy_context);
6478
CeedQFunctionContextSetData(darcy_context, CEED_MEM_HOST, CEED_COPY_VALUES,
6579
sizeof(*darcy_ctx), darcy_ctx);
6680
problem_data->qfunction_context = darcy_context;
6781
CeedQFunctionContextSetDataDestroy(darcy_context, CEED_MEM_HOST,
6882
FreeContextPetsc);
69-
83+
PetscCall( PetscFree(darcy_ctx) );
7084
PetscFunctionReturn(0);
7185
}

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

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,8 @@ PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx) {
3232
// 3) darcy3d-prism
3333

3434
// 4) richard
35-
35+
PetscCall( PetscFunctionListAdd(&app_ctx->problems, "richard2d",
36+
Hdiv_RICHARD2D) );
3637
PetscFunctionReturn(0);
3738
}
3839

0 commit comments

Comments
 (0)