1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 4 #include "../qfunctions/monitor_cfl.h" 5 #include <navierstokes.h> 6 7 typedef struct { 8 OperatorApplyContext op_monitor_ctx; 9 Vec values; 10 PetscInt num_comps; 11 PetscInt tab_level; 12 PetscBool is_header_written; 13 } *MonitorCflCtx; 14 15 static PetscErrorCode MonitorCflCtxDestroy(MonitorCflCtx *monitor_ctx) { 16 MonitorCflCtx monitor_ctx_ = *monitor_ctx; 17 18 PetscFunctionBeginUser; 19 PetscCall(OperatorApplyContextDestroy(monitor_ctx_->op_monitor_ctx)); 20 PetscCall(VecDestroy(&monitor_ctx_->values)); 21 PetscCall(PetscFree(monitor_ctx_)); 22 *monitor_ctx = NULL; 23 PetscFunctionReturn(PETSC_SUCCESS); 24 } 25 26 PetscErrorCode SetupMontiorCfl(TS ts, PetscViewerAndFormat *ctx) { 27 Honee honee; 28 Ceed ceed; 29 CeedQFunction qf_monitor = NULL; 30 CeedOperator op_monitor; 31 CeedElemRestriction elem_restr_qd, elem_restr_cfl, elem_restr_q; 32 CeedBasis basis_q; 33 CeedVector q_data; 34 CeedInt num_comp_q, q_data_size; 35 PetscInt num_comp_cfl = 1, dim, tab_level; 36 MonitorCflCtx monitor_ctx; 37 CeedQFunctionContext newt_qfctx; 38 PetscBool is_ascii; 39 40 PetscFunctionBeginUser; 41 PetscCall(PetscObjectTypeCompare((PetscObject)ctx->viewer, PETSCVIEWERASCII, &is_ascii)); 42 PetscCheck(is_ascii, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Only supports ASCII viewers"); 43 PetscCall(TSGetApplicationContext(ts, &honee)); 44 PetscCall(DMGetDimension(honee->dm, &dim)); 45 ceed = honee->ceed; 46 47 PetscCall(PetscNew(&monitor_ctx)); 48 49 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q)); 50 PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q)); 51 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_q, &num_comp_q)); 52 PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, 53 &q_data, &q_data_size)); 54 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, num_comp_cfl, &elem_restr_cfl)); 55 56 { // Get restriction and basis from the RHS function 57 CeedOperator *sub_ops, main_op = honee->op_ifunction ? honee->op_ifunction : honee->op_rhs_ctx->op; 58 PetscInt sub_op_index = 0; // will be 0 for the volume op 59 60 PetscCallCeed(ceed, CeedOperatorCompositeGetSubList(main_op, &sub_ops)); 61 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &newt_qfctx)); 62 } 63 64 switch (dim) { 65 case 2: 66 switch (honee->phys->state_var) { 67 case STATEVAR_PRIMITIVE: 68 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Prim, MonitorCFL_2D_Prim_loc, &qf_monitor)); 69 break; 70 case STATEVAR_CONSERVATIVE: 71 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Conserv, MonitorCFL_2D_Conserv_loc, &qf_monitor)); 72 break; 73 case STATEVAR_ENTROPY: 74 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Entropy, MonitorCFL_2D_Entropy_loc, &qf_monitor)); 75 break; 76 } 77 break; 78 case 3: 79 switch (honee->phys->state_var) { 80 case STATEVAR_PRIMITIVE: 81 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Prim, MonitorCFL_3D_Prim_loc, &qf_monitor)); 82 break; 83 case STATEVAR_CONSERVATIVE: 84 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Conserv, MonitorCFL_3D_Conserv_loc, &qf_monitor)); 85 break; 86 case STATEVAR_ENTROPY: 87 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Entropy, MonitorCFL_3D_Entropy_loc, &qf_monitor)); 88 break; 89 } 90 break; 91 } 92 PetscCheck(qf_monitor, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, 93 "Could not create CFL monitor QFunction for dim %" PetscInt_FMT " and state variable %s", dim, StateVariables[honee->phys->state_var]); 94 95 PetscBool is_advection; 96 PetscCall(PetscStrcmp("advection", honee->app_ctx->problem_name, &is_advection)); 97 if (is_advection) { 98 NewtonianIdealGasContext newt_ctx; 99 100 PetscCall(PetscNew(&newt_ctx)); 101 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&newt_qfctx)); 102 PetscCallCeed(ceed, CeedQFunctionContextCreate(ceed, &newt_qfctx)); 103 PetscCallCeed(ceed, CeedQFunctionContextSetData(newt_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*newt_ctx), newt_ctx)); 104 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(newt_qfctx, CEED_MEM_HOST, FreeContextPetsc)); 105 } 106 107 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_monitor, newt_qfctx)); 108 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q", num_comp_q, CEED_EVAL_INTERP)); 109 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q_data", q_data_size, CEED_EVAL_NONE)); 110 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_monitor, "v", num_comp_cfl, CEED_EVAL_NONE)); 111 112 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_monitor, NULL, NULL, &op_monitor)); 113 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 114 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data)); 115 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "v", elem_restr_cfl, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 116 117 PetscCall(OperatorApplyContextCreate(honee->dm, NULL, ceed, op_monitor, NULL, NULL, NULL, NULL, &monitor_ctx->op_monitor_ctx)); 118 119 PetscCall(CeedOperatorCreateLocalVecs(op_monitor, DMReturnVecType(honee->dm), PETSC_COMM_SELF, NULL, &monitor_ctx->values)); 120 PetscCall(VecSetBlockSize(monitor_ctx->values, num_comp_cfl)); 121 monitor_ctx->num_comps = num_comp_cfl; 122 PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level)); 123 monitor_ctx->tab_level = tab_level + 1; 124 125 ctx->data = monitor_ctx; 126 ctx->data_destroy = (PetscCtxDestroyFn *)MonitorCflCtxDestroy; 127 128 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&newt_qfctx)); 129 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 130 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 131 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_cfl)); 132 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 133 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 134 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_monitor)); 135 PetscCallCeed(ceed, CeedOperatorDestroy(&op_monitor)); 136 PetscFunctionReturn(PETSC_SUCCESS); 137 } 138 139 PetscErrorCode TSMonitor_Cfl(TS ts, PetscInt step, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) { 140 MonitorCflCtx monitor_ctx = (MonitorCflCtx)ctx->data; 141 Honee honee; 142 MPI_Comm comm; 143 TSConvergedReason reason; 144 PetscReal part_minmax[2], global_minmax[2], dt; 145 146 PetscFunctionBeginUser; 147 PetscCall(TSGetConvergedReason(ts, &reason)); 148 if (!(step % ctx->view_interval == 0 || reason != TS_CONVERGED_ITERATING)) PetscFunctionReturn(PETSC_SUCCESS); 149 PetscCall(TSGetApplicationContext(ts, &honee)); 150 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm)); 151 152 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time)); 153 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc)); 154 155 PetscCall(ApplyCeedOperatorLocalToLocal(honee->Q_loc, monitor_ctx->values, monitor_ctx->op_monitor_ctx)); 156 157 PetscCall(VecMin(monitor_ctx->values, NULL, &part_minmax[0])); 158 PetscCall(VecMax(monitor_ctx->values, NULL, &part_minmax[1])); 159 // Need to get global min/max since `values` is only the local partition 160 PetscCall(PetscGlobalMinMaxReal(comm, part_minmax, global_minmax)); 161 PetscCall(TSGetTimeStep(ts, &dt)); 162 global_minmax[0] *= dt; 163 global_minmax[1] *= dt; 164 165 if (ctx->format == PETSC_VIEWER_ASCII_CSV) { 166 if (!monitor_ctx->is_header_written) { 167 char buf[PETSC_MAX_PATH_LEN]; 168 const char *buf_const; 169 time_t t = time(NULL); 170 171 PetscCall(PetscGetVersion(buf, sizeof(buf))); 172 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_version: %s\n", buf)); 173 PetscCall(PetscGetPetscDir(&buf_const)); 174 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_directory: %s\n", buf_const)); 175 PetscCall(PetscGetArchType(buf, sizeof(buf))); 176 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_arch: %s\n", buf)); 177 PetscCheck(strftime(buf, sizeof(buf), "%FT%T%z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed"); 178 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date: %s\n", buf)); 179 PetscCheck(strftime(buf, sizeof(buf), "%Z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed"); 180 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date_timezone: %s\n", buf)); 181 PetscCall(PetscGetUserName(buf, sizeof(buf))); 182 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# username: %s\n", buf)); 183 PetscCall(PetscGetHostName(buf, sizeof(buf))); 184 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# hostname: %s\n", buf)); 185 PetscCall(PetscGetWorkingDirectory(buf, sizeof(buf))); 186 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# working_directory: %s\n", buf)); 187 188 PetscCall(PetscViewerFileGetName(ctx->viewer, &buf_const)); 189 PetscCall(PetscGetFullPath(buf_const, buf, sizeof(buf))); 190 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# original_file_path: %s\n", buf)); 191 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "#\n")); 192 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "step,time,mincfl,maxcfl\n")); 193 monitor_ctx->is_header_written = PETSC_TRUE; 194 } 195 196 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%" PetscInt_FMT ",%0.17e,", step, solution_time)); 197 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%0.17e,%0.17e", global_minmax[0], global_minmax[1])); 198 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n")); 199 } else { 200 PetscCall(PetscViewerASCIIAddTab(ctx->viewer, monitor_ctx->tab_level)); 201 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "Min, Max CFL: %0.5g, %0.5g\n", global_minmax[0], global_minmax[1])); 202 PetscCall(PetscViewerASCIISubtractTab(ctx->viewer, monitor_ctx->tab_level)); 203 } 204 PetscFunctionReturn(PETSC_SUCCESS); 205 } 206