diff --git a/examples/Hdiv-mixed/conv_test.sh b/examples/Hdiv-mixed/conv_test.sh index 1cc1488235..798f2b70eb 100755 --- a/examples/Hdiv-mixed/conv_test.sh +++ b/examples/Hdiv-mixed/conv_test.sh @@ -32,7 +32,7 @@ declare -A run_flags run_flags[pc_type]=svd if [[ $dim -eq 2 ]]; then - run_flags[problem]=darcy2d + run_flags[problem]=richard2d run_flags[dm_plex_dim]=$dim run_flags[dm_plex_box_faces]=2,2 else @@ -43,8 +43,8 @@ declare -A run_flags declare -A test_flags test_flags[res_start]=2 - test_flags[res_stride]=1 - test_flags[res_end]=12 + test_flags[res_stride]=2 + test_flags[res_end]=10 file_name=conv_test_result.csv diff --git a/examples/Hdiv-mixed/conv_test_result.csv b/examples/Hdiv-mixed/conv_test_result.csv index 9997c7c1cc..1c19c3f9ff 100644 --- a/examples/Hdiv-mixed/conv_test_result.csv +++ b/examples/Hdiv-mixed/conv_test_result.csv @@ -1,12 +1,6 @@ run,mesh_res,error_u,error_p -0,2,9.13818,0.09469 -1,3,4.60103,0.05186 -2,4,2.75883,0.03199 -3,5,1.81813,0.02132 -4,6,1.27780,0.01505 -5,7,0.93759,0.01108 -6,8,0.71073,0.00842 -7,9,0.55232,0.00655 -8,10,0.43767,0.00520 -9,11,0.35229,0.00418 -10,12,0.28723,0.00340 +0,2,8.22730,0.00852 +1,4,2.50219,0.00291 +2,6,1.60805,0.00181 +3,8,1.02228,0.00114 +4,10,0.70467,0.00078 diff --git a/examples/Hdiv-mixed/convrate_mixed.png b/examples/Hdiv-mixed/convrate_mixed.png index cd9db19a63..e98533fd9d 100644 Binary files a/examples/Hdiv-mixed/convrate_mixed.png and b/examples/Hdiv-mixed/convrate_mixed.png differ diff --git a/examples/Hdiv-mixed/include/setup-solvers.h b/examples/Hdiv-mixed/include/setup-solvers.h index a142adf840..f615d287e9 100644 --- a/examples/Hdiv-mixed/include/setup-solvers.h +++ b/examples/Hdiv-mixed/include/setup-solvers.h @@ -19,10 +19,10 @@ PetscErrorCode PDESolver(MPI_Comm comm, DM dm, Ceed ceed, CeedData ceed_data, VecType vec_type, SNES snes, KSP ksp, Vec *U); PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U, CeedScalar *l2_error_u, CeedScalar *l2_error_p); -PetscErrorCode PrintOutput(Ceed ceed, +PetscErrorCode PrintOutput(Ceed ceed, AppCtx app_ctx, PetscBool has_ts, CeedMemType mem_type_backend, - SNES snes, KSP ksp, + TS ts, SNES snes, KSP ksp, Vec U, CeedScalar l2_error_u, - CeedScalar l2_error_p, AppCtx app_ctx); + CeedScalar l2_error_p); #endif // setup_solvers_h diff --git a/examples/Hdiv-mixed/include/setup-ts.h b/examples/Hdiv-mixed/include/setup-ts.h index c8afa97e4c..4da1ea4d18 100644 --- a/examples/Hdiv-mixed/include/setup-ts.h +++ b/examples/Hdiv-mixed/include/setup-ts.h @@ -20,6 +20,6 @@ PetscErrorCode SetupResidualOperatorCtx_P0(MPI_Comm comm, DM dm, Ceed ceed, PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y, void *ctx_residual_ut); PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx, - Vec *U, PetscScalar *f_time, TS *ts); + Vec *U, TS *ts); #endif // setup_ts_h diff --git a/examples/Hdiv-mixed/include/structs.h b/examples/Hdiv-mixed/include/structs.h index c1ddd4540e..9713e485eb 100644 --- a/examples/Hdiv-mixed/include/structs.h +++ b/examples/Hdiv-mixed/include/structs.h @@ -18,8 +18,7 @@ struct AppCtx_ { // Problem type arguments PetscFunctionList problems; char problem_name[PETSC_MAX_PATH_LEN]; - CeedContextFieldLabel solution_time_label; - CeedScalar t_final; + CeedScalar t_final, t; }; // PETSc operator contexts @@ -31,8 +30,9 @@ struct OperatorApplyContext_ { CeedOperator op_apply; DM dm; Ceed ceed; - CeedScalar t; - CeedContextFieldLabel solution_time_label, final_time_label; + CeedScalar t, dt; + CeedContextFieldLabel solution_time_label, final_time_label, + timestep_label; ; }; // libCEED data struct diff --git a/examples/Hdiv-mixed/main.c b/examples/Hdiv-mixed/main.c index 1990948cae..cf8a851b5e 100644 --- a/examples/Hdiv-mixed/main.c +++ b/examples/Hdiv-mixed/main.c @@ -136,10 +136,10 @@ int main(int argc, char **argv) { //PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm, // bc_pressure) ); - // --------------------------------------------------------------------------- // Setup TSSolve for Richard problem // --------------------------------------------------------------------------- + TS ts; if (problem_data->has_ts) { // --------------------------------------------------------------------------- // Create global initial conditions @@ -154,46 +154,49 @@ int main(int argc, char **argv) { ceed_data->ctx_initial_u0, ceed_data->ctx_initial_p0, ceed_data->ctx_residual_ut); - VecView(U, PETSC_VIEWER_STDOUT_WORLD); + //VecView(U, PETSC_VIEWER_STDOUT_WORLD); + // Solve Richards problem + PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_loc) ); + PetscCall( VecZeroEntries(ceed_data->ctx_residual_ut->X_t_loc) ); + PetscCall( TSSolveRichard(dm, ceed_data, app_ctx, + &U, &ts) ); + //VecView(U, PETSC_VIEWER_STDOUT_WORLD); } + SNES snes; + KSP ksp; if (!problem_data->has_ts) { // --------------------------------------------------------------------------- - // Solve PDE + // Setup SNES for Darcy problem // --------------------------------------------------------------------------- // Create SNES - SNES snes; - KSP ksp; PetscCall( SNESCreate(comm, &snes) ); PetscCall( SNESGetKSP(snes, &ksp) ); PetscCall( PDESolver(comm, dm, ceed, ceed_data, vec_type, snes, ksp, &U) ); //VecView(U, PETSC_VIEWER_STDOUT_WORLD); + } - // --------------------------------------------------------------------------- - // Compute L2 error of mms problem - // --------------------------------------------------------------------------- - CeedScalar l2_error_u, l2_error_p; - PetscCall( ComputeL2Error(dm, ceed,ceed_data, U, &l2_error_u, - &l2_error_p) ); - - // --------------------------------------------------------------------------- - // Print output results - // --------------------------------------------------------------------------- - PetscCall( PrintOutput(ceed, mem_type_backend, - snes, ksp, U, l2_error_u, l2_error_p, app_ctx) ); - - // --------------------------------------------------------------------------- - // Save solution (paraview) - // --------------------------------------------------------------------------- - PetscViewer viewer; + // --------------------------------------------------------------------------- + // Compute L2 error of mms problem + // --------------------------------------------------------------------------- + CeedScalar l2_error_u, l2_error_p; + PetscCall( ComputeL2Error(dm, ceed, ceed_data, U, &l2_error_u, + &l2_error_p) ); - PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) ); - PetscCall( VecView(U, viewer) ); - PetscCall( PetscViewerDestroy(&viewer) ); + // --------------------------------------------------------------------------- + // Print output results + // --------------------------------------------------------------------------- + PetscCall( PrintOutput(ceed, app_ctx, problem_data->has_ts, mem_type_backend, + ts, snes, ksp, U, l2_error_u, l2_error_p) ); - PetscCall( SNESDestroy(&snes) ); + // --------------------------------------------------------------------------- + // Save solution (paraview) + // --------------------------------------------------------------------------- + PetscViewer viewer; + PetscCall( PetscViewerVTKOpen(comm,"solution.vtu",FILE_MODE_WRITE,&viewer) ); + PetscCall( VecView(U, viewer) ); + PetscCall( PetscViewerDestroy(&viewer) ); - } // --------------------------------------------------------------------------- // Free objects // --------------------------------------------------------------------------- @@ -203,6 +206,11 @@ int main(int argc, char **argv) { PetscCall( DMDestroy(&dm_u0) ); PetscCall( DMDestroy(&dm_p0) ); PetscCall( VecDestroy(&U) ); + if (problem_data->has_ts) { + PetscCall( TSDestroy(&ts) ); + } else { + PetscCall( SNESDestroy(&snes) ); + } PetscCall( CeedDataDestroy(ceed_data, problem_data) ); // -- Function list diff --git a/examples/Hdiv-mixed/problems/darcy2d.c b/examples/Hdiv-mixed/problems/darcy2d.c index c5774fd1f8..fca8b2bdd3 100644 --- a/examples/Hdiv-mixed/problems/darcy2d.c +++ b/examples/Hdiv-mixed/problems/darcy2d.c @@ -54,7 +54,7 @@ PetscErrorCode Hdiv_DARCY2D(Ceed ceed, ProblemData problem_data, void *ctx) { // ------------------------------------------------------ // Command line Options // ------------------------------------------------------ - CeedScalar kappa = 1., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.; + CeedScalar kappa = 10., rho_a0 = 998.2, g = 9.8, alpha_a = 1., b_a = 10.; PetscOptionsBegin(app_ctx->comm, NULL, "Options for Hdiv-mixed problem", NULL); PetscCall( PetscOptionsScalar("-kappa", "Hydraulic Conductivity", NULL, kappa, &kappa, NULL)); diff --git a/examples/Hdiv-mixed/problems/richard2d.c b/examples/Hdiv-mixed/problems/richard2d.c index 0df463f217..7240a5375a 100644 --- a/examples/Hdiv-mixed/problems/richard2d.c +++ b/examples/Hdiv-mixed/problems/richard2d.c @@ -21,6 +21,7 @@ #include "../qfunctions/richard-system2d.h" #include "../qfunctions/richard-true2d.h" #include "../qfunctions/richard-ics2d.h" +#include "../qfunctions/darcy-error2d.h" //#include "../qfunctions/pressure-boundary2d.h" #include "petscsystypes.h" @@ -50,14 +51,12 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { problem_data->rhs_p0_loc = RichardRhsP02D_loc; problem_data->ics_p = RichardICsP2D; problem_data->ics_p_loc = RichardICsP2D_loc; - //problem_data->ics_p = RichardICsP2D; - //problem_data->ics_p_loc = RichardICsP2D_loc; - //problem_data->residual = RichardSystem2D; - //problem_data->residual_loc = RichardSystem2D_loc; + problem_data->residual = RichardSystem2D; + problem_data->residual_loc = RichardSystem2D_loc; //problem_data->jacobian = JacobianRichardSystem2D; //problem_data->jacobian_loc = JacobianRichardSystem2D_loc; - //problem_data->error = DarcyError2D; - //problem_data->error_loc = DarcyError2D_loc; + problem_data->error = DarcyError2D; + problem_data->error_loc = DarcyError2D_loc; //problem_data->bc_pressure = BCPressure2D; //problem_data->bc_pressure_loc = BCPressure2D_loc; problem_data->has_ts = PETSC_TRUE; @@ -81,7 +80,7 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { rho_a0, &rho_a0, NULL)); PetscCall( PetscOptionsScalar("-beta", "Water compressibility", NULL, beta, &beta, NULL)); - app_ctx->t_final = 5.; + app_ctx->t_final = 0.5; PetscCall( PetscOptionsScalar("-t_final", "End time", NULL, app_ctx->t_final, &app_ctx->t_final, NULL)); PetscOptionsEnd(); @@ -105,10 +104,15 @@ PetscErrorCode Hdiv_RICHARD2D(Ceed ceed, ProblemData problem_data, void *ctx) { offsetof(struct RICHARDContext_, t), 1, "current solver time"); CeedQFunctionContextRegisterDouble(richard_context, "final_time", offsetof(struct RICHARDContext_, t_final), 1, "final time"); + CeedQFunctionContextRegisterDouble(richard_context, "time_step", + offsetof(struct RICHARDContext_, dt), 1, "time step"); problem_data->true_qfunction_ctx = richard_context; CeedQFunctionContextReferenceCopy(richard_context, &problem_data->rhs_u0_qfunction_ctx); - + CeedQFunctionContextReferenceCopy(richard_context, + &problem_data->residual_qfunction_ctx); + CeedQFunctionContextReferenceCopy(richard_context, + &problem_data->error_qfunction_ctx); PetscCall( PetscFree(richard_ctx) ); PetscFunctionReturn(0); diff --git a/examples/Hdiv-mixed/qfunctions/darcy-system2d.h b/examples/Hdiv-mixed/qfunctions/darcy-system2d.h index 98df9aa3ff..07edd820e9 100644 --- a/examples/Hdiv-mixed/qfunctions/darcy-system2d.h +++ b/examples/Hdiv-mixed/qfunctions/darcy-system2d.h @@ -172,7 +172,7 @@ CEED_QFUNCTION(JacobianDarcySystem2D)(void *ctx, CeedInt Q, // *INDENT-ON* DARCYContext context = (DARCYContext)ctx; - const CeedScalar kappa = context->kappa; + const CeedScalar kappa = context->kappa; const CeedScalar rho_a0 = context->rho_a0; const CeedScalar g = context->g; const CeedScalar alpha_a = context->alpha_a; diff --git a/examples/Hdiv-mixed/qfunctions/darcy-true2d.h b/examples/Hdiv-mixed/qfunctions/darcy-true2d.h index 5a85f50559..3eaef40f0f 100644 --- a/examples/Hdiv-mixed/qfunctions/darcy-true2d.h +++ b/examples/Hdiv-mixed/qfunctions/darcy-true2d.h @@ -77,17 +77,20 @@ CEED_QFUNCTION(DarcyTrue2D)(void *ctx, const CeedInt Q, CeedScalar x = coords[i+0*Q], y = coords[i+1*Q]; CeedScalar psi = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y); CeedScalar psi_x = PI_DOUBLE*cos(PI_DOUBLE*x)*sin(PI_DOUBLE*y); + CeedScalar psi_xx = -PI_DOUBLE*PI_DOUBLE*psi; CeedScalar psi_y = PI_DOUBLE*sin(PI_DOUBLE*x)*cos(PI_DOUBLE*y); - + CeedScalar psi_yy = -PI_DOUBLE*PI_DOUBLE*psi; // k_r = b_a + alpha_a * (1 - x*y) CeedScalar k_r = b_a + alpha_a*(1-x*y); + CeedScalar k_rx = -alpha_a*y; + CeedScalar k_ry = -alpha_a*x; // rho = rho_a/rho_a0 CeedScalar rho = 1.; // u = -rho*k_r*K *[grad(\psi) - rho*g_u] - CeedScalar u[2] = {-rho*k_r*kappa*psi_x, -rho*k_r*kappa*(psi_y-1)}; - CeedScalar div_u = -rho*kappa*(-alpha_a*y*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi - -alpha_a*x*psi_y - k_r*PI_DOUBLE*PI_DOUBLE*psi); - + CeedScalar u[2] = {-rho*kappa*k_r*psi_x, + -rho*kappa*k_r*(psi_y-1)}; + CeedScalar div_u = -rho*kappa*(k_rx*psi_x + k_r*psi_xx + + k_ry*(psi_y-1) + k_r*psi_yy); // True Force: f = \div(u) true_force[i+0*Q] = div_u; // True Solution diff --git a/examples/Hdiv-mixed/qfunctions/darcy-true3d.h b/examples/Hdiv-mixed/qfunctions/darcy-true3d.h index 9c76857886..dd18bdd998 100644 --- a/examples/Hdiv-mixed/qfunctions/darcy-true3d.h +++ b/examples/Hdiv-mixed/qfunctions/darcy-true3d.h @@ -77,20 +77,26 @@ CEED_QFUNCTION(DarcyTrue3D)(void *ctx, const CeedInt Q, CeedScalar x = coords[i+0*Q], y = coords[i+1*Q], z = coords[i+2*Q]; CeedScalar psi = sin(PI_DOUBLE*x) * sin(PI_DOUBLE*y) * sin(PI_DOUBLE*z); CeedScalar psi_x = PI_DOUBLE*cos(PI_DOUBLE*x) *sin(PI_DOUBLE*y) *sin(PI_DOUBLE*z); + CeedScalar psi_xx = -PI_DOUBLE*PI_DOUBLE*psi; CeedScalar psi_y = PI_DOUBLE*sin(PI_DOUBLE*x) *cos(PI_DOUBLE*y) *sin(PI_DOUBLE*z); + CeedScalar psi_yy = -PI_DOUBLE*PI_DOUBLE*psi; CeedScalar psi_z = PI_DOUBLE*sin(PI_DOUBLE*x) *sin(PI_DOUBLE*y) *cos(PI_DOUBLE*z); + CeedScalar psi_zz = -PI_DOUBLE*PI_DOUBLE*psi; // k_r = b_a + alpha_a * (psi - x2) CeedScalar k_r = b_a + alpha_a*(1-x*y*z); + CeedScalar k_rx = -alpha_a*y*z; + CeedScalar k_ry = -alpha_a*x*z; + CeedScalar k_rz = -alpha_a*x*y; // rho = rho_a/rho_a0 CeedScalar rho = 1.; // u = -rho*k_r*K *[grad(\psi) - rho*g_u] - CeedScalar u[3] = {-rho*k_r*kappa*psi_x, - -rho*k_r*kappa*psi_y, - -rho*k_r*kappa*(psi_z-1)}; - CeedScalar div_u = -rho*kappa*(-alpha_a*y*z*psi_x - k_r*PI_DOUBLE*PI_DOUBLE*psi - -alpha_a*x*z*psi_y - k_r*PI_DOUBLE*PI_DOUBLE*psi - -alpha_a*x*y*psi_z - k_r*PI_DOUBLE*PI_DOUBLE*psi); + CeedScalar u[3] = {-rho*kappa*k_r*psi_x, + -rho*kappa*k_r*psi_y, + -rho*kappa*k_r*(psi_z-1)}; + CeedScalar div_u = -rho*kappa*(k_rx*psi_x + k_r*psi_xx + + k_ry*psi_y + k_r*psi_yy + + k_rz*(psi_z-1) + k_r*psi_zz); // True Force: f = \div(u) true_force[i+0*Q] = div_u; diff --git a/examples/Hdiv-mixed/qfunctions/richard-ics2d.h b/examples/Hdiv-mixed/qfunctions/richard-ics2d.h index a50bce6663..49664a40ae 100644 --- a/examples/Hdiv-mixed/qfunctions/richard-ics2d.h +++ b/examples/Hdiv-mixed/qfunctions/richard-ics2d.h @@ -66,7 +66,7 @@ struct RICHARDContext_ { CeedScalar rho_a0; CeedScalar alpha_a, b_a; CeedScalar beta, p0; - CeedScalar t, t_final; + CeedScalar t, t_final, dt; CeedScalar gamma; }; #endif @@ -105,8 +105,9 @@ CEED_QFUNCTION(RichardRhsU02D)(void *ctx, const CeedInt Q, CeedScalar k_r = b_a + alpha_a*(1-x*y); // rho = rho_a/rho_a0 CeedScalar rho = 1.; - // ue = -rho*k_r*K *[grad(\psi) - rho*g_u] - CeedScalar ue[2] = {-rho*k_r*kappa*psi1_x, -rho*k_r*kappa*(psi1_y-1)}; + // ue = -rho*k_r*K *[grad(\psi)] + CeedScalar ue[2] = {-rho*k_r*kappa*psi1_x, + -rho*k_r*kappa*psi1_y}; CeedScalar rhs1[2]; // rhs = (v, ue) = J^T*ue*w AlphaMatTransposeVecMult2x2(w[i], J, ue, rhs1); @@ -184,10 +185,10 @@ CEED_QFUNCTION(RichardRhsP02D)(void *ctx, const CeedInt Q, {dxdX[0][1][i], dxdX[1][1][i]}}; const CeedScalar det_J = MatDet2x2(J); // psi = exp(-gamma*t)*sin(pi*x)*sin(pi*y) - CeedScalar psi_e = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y); + CeedScalar psi1 = sin(PI_DOUBLE*x)*sin(PI_DOUBLE*y); // rhs = (q, pe) = pe*w*det_J - rhs_p0[i] = psi_e*w[i]*det_J; + rhs_p0[i] = psi1*w[i]*det_J; } // End of Quadrature Point Loop return 0; } diff --git a/examples/Hdiv-mixed/qfunctions/richard-system2d.h b/examples/Hdiv-mixed/qfunctions/richard-system2d.h index 1540b9a6ff..8dc9f63bc2 100644 --- a/examples/Hdiv-mixed/qfunctions/richard-system2d.h +++ b/examples/Hdiv-mixed/qfunctions/richard-system2d.h @@ -22,6 +22,7 @@ #include #include +#include "ceed/ceed-f64.h" #include "utils.h" // See Matthew Farthing, Christopher Kees, Cass Miller (2003) @@ -72,14 +73,13 @@ struct RICHARDContext_ { CeedScalar rho_a0; CeedScalar alpha_a, b_a; CeedScalar beta, p0; - CeedScalar t, t_final; + CeedScalar t, t_final, dt; CeedScalar gamma; }; #endif // ----------------------------------------------------------------------------- // Residual evaluation for Richard problem // ----------------------------------------------------------------------------- -/* CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { @@ -90,8 +90,9 @@ CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, (*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[2], (*div_u) = (const CeedScalar(*))in[3], (*p) = (const CeedScalar(*))in[4], - (*p_t) = (const CeedScalar(*))in[5], - (*coords) = in[6]; + (*f) = in[5], + (*coords) = in[6], + (*p_t) = (const CeedScalar(*))in[7]; // Outputs CeedScalar (*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0], @@ -99,16 +100,18 @@ CEED_QFUNCTION(RichardSystem2D)(void *ctx, CeedInt Q, (*q) = (CeedScalar(*))out[2]; // Context RICHARDContext context = (RICHARDContext)ctx; - const CeedScalar kappa = context->kappa; - const CeedScalar alpha_a = context->alpha_a; - const CeedScalar b_a = context->b_a; + const CeedScalar kappa = context->kappa; const CeedScalar rho_a0 = context->rho_a0; - const CeedScalar beta = context->beta; - const CeedScalar g = context->g; - const CeedScalar p0 = context->p0; // atmospheric pressure - CeedScalar t = context->time; - // *INDENT-ON* + const CeedScalar g = context->g; + const CeedScalar alpha_a = context->alpha_a; + const CeedScalar b_a = context->b_a; + //const CeedScalar beta = context->beta; + //const CeedScalar p0 = context->p0; // atmospheric pressure + const CeedScalar gamma = context->gamma; + CeedScalar t = context->t; + //CeedScalar dt = context->dt; + // *INDENT-ON* // Quadrature Point Loop CeedPragmaSIMD for (CeedInt i=0; ihas_ts) { - CeedVectorDestroy(&ceed_data->x_ceed); - CeedVectorDestroy(&ceed_data->y_ceed); - } + CeedVectorDestroy(&ceed_data->x_ceed); + CeedVectorDestroy(&ceed_data->y_ceed); + CeedVectorDestroy(&ceed_data->x_t_ceed); CeedVectorDestroy(&ceed_data->x_coord); // Restrictions CeedElemRestrictionDestroy(&ceed_data->elem_restr_x); @@ -55,16 +54,17 @@ PetscErrorCode CeedDataDestroy(CeedData ceed_data, ProblemData problem_data) { CeedQFunctionDestroy(&ceed_data->qf_ics_p); CeedOperatorDestroy(&ceed_data->op_ics_p); } - + //QFunctions + CeedQFunctionDestroy(&ceed_data->qf_residual); + CeedQFunctionDestroy(&ceed_data->qf_error); + //Operators + CeedOperatorDestroy(&ceed_data->op_residual); + CeedOperatorDestroy(&ceed_data->op_error); if (!problem_data->has_ts) { //QFunctions - CeedQFunctionDestroy(&ceed_data->qf_residual); CeedQFunctionDestroy(&ceed_data->qf_jacobian); - CeedQFunctionDestroy(&ceed_data->qf_error); //Operators - CeedOperatorDestroy(&ceed_data->op_residual); CeedOperatorDestroy(&ceed_data->op_jacobian); - CeedOperatorDestroy(&ceed_data->op_error); } PetscCall( PetscFree(ceed_data) ); @@ -485,63 +485,80 @@ PetscErrorCode SetupLibceed(DM dm, DM dm_u0, DM dm_p0, Ceed ceed, NULL); } - if (!problem_data->has_ts) { - // --------------------------------------------------------------------------- - // Persistent libCEED vectors - // --------------------------------------------------------------------------- - // -- Operator action variables: we use them in setup-solvers.c - CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->x_ceed, - NULL); - CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->y_ceed, - NULL); - - // Local residual evaluator - // --------------------------------------------------------------------------- - // Create the QFunction and Operator that computes the residual of the PDE. - // --------------------------------------------------------------------------- - // -- QFunction - CeedQFunctionCreateInterior(ceed, 1, problem_data->residual, - problem_data->residual_loc, &qf_residual); - CeedQFunctionSetContext(qf_residual, problem_data->residual_qfunction_ctx); - CeedQFunctionContextDestroy(&problem_data->residual_qfunction_ctx); - CeedQFunctionAddInput(qf_residual, "weight", 1, CEED_EVAL_WEIGHT); - CeedQFunctionAddInput(qf_residual, "dx", dim*dim, CEED_EVAL_GRAD); - CeedQFunctionAddInput(qf_residual, "u", dim, CEED_EVAL_INTERP); - CeedQFunctionAddInput(qf_residual, "div_u", 1, CEED_EVAL_DIV); - CeedQFunctionAddInput(qf_residual, "p", 1, CEED_EVAL_INTERP); - CeedQFunctionAddInput(qf_residual, "true force", 1, CEED_EVAL_NONE); - CeedQFunctionAddInput(qf_residual, "x", num_comp_x, CEED_EVAL_INTERP); - CeedQFunctionAddOutput(qf_residual, "v", dim, CEED_EVAL_INTERP); - CeedQFunctionAddOutput(qf_residual, "div_v", 1, CEED_EVAL_DIV); - CeedQFunctionAddOutput(qf_residual, "q", 1, CEED_EVAL_INTERP); - - // -- Operator - CeedOperatorCreate(ceed, qf_residual, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, - &op_residual); - CeedOperatorSetField(op_residual, "weight", CEED_ELEMRESTRICTION_NONE, - ceed_data->basis_x, CEED_VECTOR_NONE); - CeedOperatorSetField(op_residual, "dx", ceed_data->elem_restr_x, - ceed_data->basis_x, ceed_data->x_coord); - CeedOperatorSetField(op_residual, "u", ceed_data->elem_restr_u, - ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_residual, "div_u", ceed_data->elem_restr_u, - ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_residual, "p", ceed_data->elem_restr_p, - ceed_data->basis_p, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_residual, "true force", ceed_data->elem_restr_p_i, - CEED_BASIS_COLLOCATED, true_force); - CeedOperatorSetField(op_residual, "x", ceed_data->elem_restr_x, - ceed_data->basis_x, ceed_data->x_coord); - CeedOperatorSetField(op_residual, "v", ceed_data->elem_restr_u, - ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_residual, "div_v", ceed_data->elem_restr_u, - ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_residual, "q", ceed_data->elem_restr_p, - ceed_data->basis_p, CEED_VECTOR_ACTIVE); - // -- Save libCEED data to apply operator in matops.c - ceed_data->qf_residual = qf_residual; - ceed_data->op_residual = op_residual; + // --------------------------------------------------------------------------- + // Persistent libCEED vectors + // --------------------------------------------------------------------------- + // -- Operator action variables: we use them in setup-solvers.c/setup-ts.c + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->x_ceed, + NULL); + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->y_ceed, + NULL); + // -- Operator action variables: we use them in setup-ts.c + CeedElemRestrictionCreateVector(ceed_data->elem_restr_u, &ceed_data->x_t_ceed, + NULL); + // Local residual evaluator + // --------------------------------------------------------------------------- + // Create the QFunction and Operator that computes the residual of the PDE. + // --------------------------------------------------------------------------- + // -- QFunction + CeedQFunctionCreateInterior(ceed, 1, problem_data->residual, + problem_data->residual_loc, &qf_residual); + CeedQFunctionSetContext(qf_residual, problem_data->residual_qfunction_ctx); + CeedQFunctionContextDestroy(&problem_data->residual_qfunction_ctx); + CeedQFunctionAddInput(qf_residual, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddInput(qf_residual, "dx", dim*dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_residual, "u", dim, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_residual, "div_u", 1, CEED_EVAL_DIV); + CeedQFunctionAddInput(qf_residual, "p", 1, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_residual, "true force", 1, CEED_EVAL_NONE); + CeedQFunctionAddInput(qf_residual, "x", num_comp_x, CEED_EVAL_INTERP); + if (problem_data->has_ts) { + CeedQFunctionAddInput(qf_residual, "p_t", 1, CEED_EVAL_INTERP); + } + CeedQFunctionAddOutput(qf_residual, "v", dim, CEED_EVAL_INTERP); + CeedQFunctionAddOutput(qf_residual, "div_v", 1, CEED_EVAL_DIV); + CeedQFunctionAddOutput(qf_residual, "q", 1, CEED_EVAL_INTERP); + // -- Operator + CeedOperatorCreate(ceed, qf_residual, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, + &op_residual); + if (problem_data->has_ts) { + //double t = ceed_data->ctx_residual_ut->t; + CeedOperatorContextGetFieldLabel(op_residual, "time", + &ceed_data->ctx_residual_ut->solution_time_label); + //CeedOperatorContextGetFieldLabel(op_residual, "time_step", + // &ceed_data->ctx_residual_ut->timestep_label); + //CeedOperatorContextSetDouble(op_residual, + // ceed_data->ctx_residual_ut->solution_time_label, &t); + } + CeedOperatorSetField(op_residual, "weight", CEED_ELEMRESTRICTION_NONE, + ceed_data->basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_residual, "dx", ceed_data->elem_restr_x, + ceed_data->basis_x, ceed_data->x_coord); + CeedOperatorSetField(op_residual, "u", ceed_data->elem_restr_u, + ceed_data->basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_residual, "div_u", ceed_data->elem_restr_u, + ceed_data->basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_residual, "p", ceed_data->elem_restr_p, + ceed_data->basis_p, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_residual, "true force", ceed_data->elem_restr_p_i, + CEED_BASIS_COLLOCATED, true_force); + CeedOperatorSetField(op_residual, "x", ceed_data->elem_restr_x, + ceed_data->basis_x, ceed_data->x_coord); + if (problem_data->has_ts) { + CeedOperatorSetField(op_residual, "p_t", ceed_data->elem_restr_p, + ceed_data->basis_p, ceed_data->x_t_ceed); + } + CeedOperatorSetField(op_residual, "v", ceed_data->elem_restr_u, + ceed_data->basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_residual, "div_v", ceed_data->elem_restr_u, + ceed_data->basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_residual, "q", ceed_data->elem_restr_p, + ceed_data->basis_p, CEED_VECTOR_ACTIVE); + // -- Save libCEED data to apply operator in matops.c + ceed_data->qf_residual = qf_residual; + ceed_data->op_residual = op_residual; + if (!problem_data->has_ts) { // --------------------------------------------------------------------------- // Add Pressure boundary condition. See setup-boundary.c // --------------------------------------------------------------------------- @@ -595,40 +612,40 @@ PetscErrorCode SetupLibceed(DM dm, DM dm_u0, DM dm_p0, Ceed ceed, // -- Save libCEED data to apply operator in matops.c ceed_data->qf_jacobian = qf_jacobian; ceed_data->op_jacobian = op_jacobian; - - // --------------------------------------------------------------------------- - // Setup Error Qfunction - // --------------------------------------------------------------------------- - // Create the q-function that sets up the error - CeedQFunctionCreateInterior(ceed, 1, problem_data->error, - problem_data->error_loc, &qf_error); - CeedQFunctionSetContext(qf_error, problem_data->error_qfunction_ctx); - CeedQFunctionContextDestroy(&problem_data->error_qfunction_ctx); - CeedQFunctionAddInput(qf_error, "weight", 1, CEED_EVAL_WEIGHT); - CeedQFunctionAddInput(qf_error, "dx", dim*dim, CEED_EVAL_GRAD); - CeedQFunctionAddInput(qf_error, "u", dim, CEED_EVAL_INTERP); - CeedQFunctionAddInput(qf_error, "p", 1, CEED_EVAL_INTERP); - CeedQFunctionAddInput(qf_error, "true solution", dim+1, CEED_EVAL_NONE); - CeedQFunctionAddOutput(qf_error, "error", dim+1, CEED_EVAL_NONE); - // Create the operator that builds the error - CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, - &op_error); - CeedOperatorSetField(op_error, "weight", CEED_ELEMRESTRICTION_NONE, - ceed_data->basis_x, CEED_VECTOR_NONE); - CeedOperatorSetField(op_error, "dx", ceed_data->elem_restr_x, - ceed_data->basis_x, ceed_data->x_coord); - CeedOperatorSetField(op_error, "u", ceed_data->elem_restr_u, - ceed_data->basis_u, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_error, "p", ceed_data->elem_restr_p, - ceed_data->basis_p, CEED_VECTOR_ACTIVE); - CeedOperatorSetField(op_error, "true solution", ceed_data->elem_restr_U_i, - CEED_BASIS_COLLOCATED, true_vec); - CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_U_i, - CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); - // -- Save libCEED data to apply operator in matops.c - ceed_data->qf_error = qf_error; - ceed_data->op_error = op_error; } + // --------------------------------------------------------------------------- + // Setup Error Qfunction + // --------------------------------------------------------------------------- + // Create the q-function that sets up the error + CeedQFunctionCreateInterior(ceed, 1, problem_data->error, + problem_data->error_loc, &qf_error); + CeedQFunctionSetContext(qf_error, problem_data->error_qfunction_ctx); + CeedQFunctionContextDestroy(&problem_data->error_qfunction_ctx); + CeedQFunctionAddInput(qf_error, "weight", 1, CEED_EVAL_WEIGHT); + CeedQFunctionAddInput(qf_error, "dx", dim*dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_error, "u", dim, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_error, "p", 1, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_error, "true solution", dim+1, CEED_EVAL_NONE); + CeedQFunctionAddOutput(qf_error, "error", dim+1, CEED_EVAL_NONE); + // Create the operator that builds the error + CeedOperatorCreate(ceed, qf_error, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, + &op_error); + CeedOperatorSetField(op_error, "weight", CEED_ELEMRESTRICTION_NONE, + ceed_data->basis_x, CEED_VECTOR_NONE); + CeedOperatorSetField(op_error, "dx", ceed_data->elem_restr_x, + ceed_data->basis_x, ceed_data->x_coord); + CeedOperatorSetField(op_error, "u", ceed_data->elem_restr_u, + ceed_data->basis_u, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_error, "p", ceed_data->elem_restr_p, + ceed_data->basis_p, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_error, "true solution", ceed_data->elem_restr_U_i, + CEED_BASIS_COLLOCATED, true_vec); + CeedOperatorSetField(op_error, "error", ceed_data->elem_restr_U_i, + CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); + // -- Save libCEED data to apply operator in matops.c + ceed_data->qf_error = qf_error; + ceed_data->op_error = op_error; + // -- Cleanup CeedVectorDestroy(&true_vec); CeedVectorDestroy(&true_force); diff --git a/examples/Hdiv-mixed/src/setup-solvers.c b/examples/Hdiv-mixed/src/setup-solvers.c index 8defce4d69..4f47c4fcc4 100644 --- a/examples/Hdiv-mixed/src/setup-solvers.c +++ b/examples/Hdiv-mixed/src/setup-solvers.c @@ -264,11 +264,10 @@ PetscErrorCode ComputeL2Error(DM dm, Ceed ceed, CeedData ceed_data, Vec U, // ----------------------------------------------------------------------------- // This function print the output // ----------------------------------------------------------------------------- -PetscErrorCode PrintOutput(Ceed ceed, +PetscErrorCode PrintOutput(Ceed ceed, AppCtx app_ctx, PetscBool has_ts, CeedMemType mem_type_backend, - SNES snes, KSP ksp, - Vec U, CeedScalar l2_error_u, - CeedScalar l2_error_p, AppCtx app_ctx) { + TS ts, SNES snes, KSP ksp, + Vec U, CeedScalar l2_error_u, CeedScalar l2_error_p) { PetscFunctionBeginUser; @@ -305,6 +304,25 @@ PetscErrorCode PrintOutput(Ceed ceed, " Owned nodes (u + p) : %" PetscInt_FMT "\n", app_ctx->problem_name, U_g_size, U_l_size ) ); + // --TS + if (has_ts) { + PetscInt ts_steps; + TSType ts_type; + TSConvergedReason ts_reason; + PetscCall( TSGetStepNumber(ts, &ts_steps) ); + PetscCall( TSGetType(ts, &ts_type) ); + PetscCall( TSGetConvergedReason(ts, &ts_reason) ); + PetscCall( PetscPrintf(app_ctx->comm, + " TS:\n" + " TS Type : %s\n" + " TS Convergence : %s\n" + " Number of TS steps : %" PetscInt_FMT "\n" + " Final time : %g\n", + ts_type, TSConvergedReasons[ts_reason], + ts_steps, (double)app_ctx->t_final) ); + + PetscCall( TSGetSNES(ts, &snes) ); + } // -- SNES PetscInt its, snes_its = 0; PetscCall( SNESGetIterationNumber(snes, &its) ); @@ -323,30 +341,31 @@ PetscErrorCode PrintOutput(Ceed ceed, " Final rnorm : %e\n", snes_type, SNESConvergedReasons[snes_reason], snes_its, (double)snes_rnorm) ); - - PetscInt ksp_its = 0; - PetscCall( SNESGetLinearSolveIterations(snes, &its) ); - ksp_its += its; - KSPType ksp_type; - KSPConvergedReason ksp_reason; - PetscReal ksp_rnorm; - PC pc; - PCType pc_type; - PetscCall( KSPGetPC(ksp, &pc) ); - PetscCall( PCGetType(pc, &pc_type) ); - PetscCall( KSPGetType(ksp, &ksp_type) ); - PetscCall( KSPGetConvergedReason(ksp, &ksp_reason) ); - PetscCall( KSPGetIterationNumber(ksp, &ksp_its) ); - PetscCall( KSPGetResidualNorm(ksp, &ksp_rnorm) ); - PetscCall( PetscPrintf(app_ctx->comm, - " KSP:\n" - " KSP Type : %s\n" - " PC Type : %s\n" - " KSP Convergence : %s\n" - " Total KSP Iterations : %" PetscInt_FMT "\n" - " Final rnorm : %e\n", - ksp_type, pc_type, KSPConvergedReasons[ksp_reason], ksp_its, - (double)ksp_rnorm ) ); + if (!has_ts) { + PetscInt ksp_its = 0; + PetscCall( SNESGetLinearSolveIterations(snes, &its) ); + ksp_its += its; + KSPType ksp_type; + KSPConvergedReason ksp_reason; + PetscReal ksp_rnorm; + PC pc; + PCType pc_type; + PetscCall( KSPGetPC(ksp, &pc) ); + PetscCall( PCGetType(pc, &pc_type) ); + PetscCall( KSPGetType(ksp, &ksp_type) ); + PetscCall( KSPGetConvergedReason(ksp, &ksp_reason) ); + PetscCall( KSPGetIterationNumber(ksp, &ksp_its) ); + PetscCall( KSPGetResidualNorm(ksp, &ksp_rnorm) ); + PetscCall( PetscPrintf(app_ctx->comm, + " KSP:\n" + " KSP Type : %s\n" + " PC Type : %s\n" + " KSP Convergence : %s\n" + " Total KSP Iterations : %" PetscInt_FMT "\n" + " Final rnorm : %e\n", + ksp_type, pc_type, KSPConvergedReasons[ksp_reason], ksp_its, + (double)ksp_rnorm ) ); + } PetscCall( PetscPrintf(app_ctx->comm, " L2 Error (MMS):\n" diff --git a/examples/Hdiv-mixed/src/setup-ts.c b/examples/Hdiv-mixed/src/setup-ts.c index db6afb086a..dd0728951b 100644 --- a/examples/Hdiv-mixed/src/setup-ts.c +++ b/examples/Hdiv-mixed/src/setup-ts.c @@ -250,8 +250,18 @@ PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y, PetscFunctionBeginUser; // Update time dependent data - - + if(ctx->t != time) { + CeedOperatorContextSetDouble(ctx->op_apply, + ctx->solution_time_label, &time); + ctx->t = time; + } + //PetscScalar dt; + //PetscCall( TSGetTimeStep(ts, &dt) ); + //if (ctx->dt != dt) { + // CeedOperatorContextSetDouble(ctx->op_apply, + // ctx->timestep_label, &dt); + // ctx->dt = dt; + //} // Global-to-local PetscCall( DMGlobalToLocal(ctx->dm, X, INSERT_VALUES, ctx->X_loc) ); PetscCall( DMGlobalToLocal(ctx->dm, X_t, INSERT_VALUES, ctx->X_t_loc) ); @@ -282,31 +292,30 @@ PetscErrorCode TSFormIResidual(TS ts, PetscReal time, Vec X, Vec X_t, Vec Y, PetscCall( VecZeroEntries(Y) ); PetscCall( DMLocalToGlobal(ctx->dm, ctx->Y_loc, ADD_VALUES, Y) ); - // Restore vectors - PetscCall( DMRestoreLocalVector(ctx->dm, &ctx->Y_loc) ); - PetscFunctionReturn(0); } // TS: Create, setup, and solve PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx, - Vec *U, PetscScalar *f_time, TS *ts) { + Vec *U, TS *ts) { MPI_Comm comm = app_ctx->comm; TSAdapt adapt; PetscFunctionBeginUser; PetscCall( TSCreate(comm, ts) ); PetscCall( TSSetDM(*ts, dm) ); + PetscCall( TSSetType(*ts, TSBDF) ); PetscCall( TSSetIFunction(*ts, NULL, TSFormIResidual, ceed_data->ctx_residual_ut) ); - PetscCall( TSSetMaxTime(*ts, 10) ); + PetscCall( TSSetMaxTime(*ts, app_ctx->t_final) ); PetscCall( TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER) ); PetscCall( TSSetTimeStep(*ts, 1.e-2) ); PetscCall( TSGetAdapt(*ts, &adapt) ); PetscCall( TSAdaptSetStepLimits(adapt, 1.e-12, 1.e2) ); PetscCall( TSSetFromOptions(*ts) ); - + ceed_data->ctx_residual_ut->t = -1.0; + //ceed_data->ctx_residual_ut->dt = -1.0; // Solve PetscScalar start_time; PetscCall( TSGetTime(*ts, &start_time) ); @@ -319,7 +328,7 @@ PetscErrorCode TSSolveRichard(DM dm, CeedData ceed_data, AppCtx app_ctx, PetscScalar final_time; PetscCall( TSGetSolveTime(*ts, &final_time) ); - *f_time = final_time; + app_ctx->t_final = final_time; PetscFunctionReturn(0); }