1 // SPDX-FileCopyrightText: Copyright (c) 2017-2024, HONEE contributors. 2 // SPDX-License-Identifier: Apache-2.0 OR BSD-2-Clause 3 #include "../qfunctions/monitor_totalkineticenergy.h" 4 #include <honee.h> 5 #include <navierstokes.h> 6 #include <time.h> 7 8 typedef struct { 9 OperatorApplyContext op_monitor_ctx; 10 Vec values; 11 PetscScalar *sum_values; 12 PetscScalar volume; 13 PetscInt num_comps; 14 PetscInt tab_level; 15 PetscBool is_header_written; 16 } *MonitorTotalKE; 17 18 static PetscErrorCode MonitorTotalKEDestroy(MonitorTotalKE *monitor_ctx) { 19 MonitorTotalKE monitor_ctx_ = *monitor_ctx; 20 21 PetscFunctionBeginUser; 22 PetscCall(OperatorApplyContextDestroy(monitor_ctx_->op_monitor_ctx)); 23 PetscCall(VecDestroy(&monitor_ctx_->values)); 24 PetscCall(PetscFree(monitor_ctx_->sum_values)); 25 PetscCall(PetscFree(monitor_ctx_)); 26 PetscFunctionReturn(PETSC_SUCCESS); 27 } 28 29 /** @brief Create CeedElemRestriction for collocated data in component-major order. 30 a. Sets the strides of the restriction to component-major order 31 Number of quadrature points is used from the CeedBasis, and number of elements is used from the CeedElemRestriction. 32 */ 33 static PetscErrorCode CreateElemRestrColloc_CompMajor(Ceed ceed, CeedInt num_comp, CeedBasis basis, CeedElemRestriction elem_restr_base, 34 CeedElemRestriction *elem_restr_collocated) { 35 CeedInt num_elem_qpts, loc_num_elem; 36 37 PetscFunctionBeginUser; 38 PetscCallCeed(ceed, CeedBasisGetNumQuadraturePoints(basis, &num_elem_qpts)); 39 PetscCallCeed(ceed, CeedElemRestrictionGetNumElements(elem_restr_base, &loc_num_elem)); 40 41 const CeedInt strides[] = {num_comp, 1, num_elem_qpts * num_comp}; 42 PetscCallCeed(ceed, CeedElemRestrictionCreateStrided(ceed, loc_num_elem, num_elem_qpts, num_comp, num_comp * loc_num_elem * num_elem_qpts, strides, 43 elem_restr_collocated)); 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 PetscErrorCode SetupMontiorTotalKineticEnergy(TS ts, PetscViewerAndFormat *ctx) { 48 Honee honee; 49 Ceed ceed; 50 CeedQFunction qf_monitor; 51 CeedOperator op_monitor; 52 CeedElemRestriction elem_restr_qd, elem_restr_totalke, elem_restr_q; 53 CeedBasis basis_q; 54 CeedVector q_data; 55 CeedInt num_comp_q, q_data_size; 56 DMLabel domain_label = NULL; 57 PetscInt label_value = 0, num_comp_totalke = 5, dim, tab_level; 58 MonitorTotalKE monitor_ctx; 59 CeedQFunctionContext newt_qfctx; 60 PetscBool is_ascii; 61 62 PetscFunctionBeginUser; 63 PetscCall(PetscObjectTypeCompare((PetscObject)ctx->viewer, PETSCVIEWERASCII, &is_ascii)); 64 PetscCheck(is_ascii, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Only supports ASCII viewers"); 65 PetscCall(TSGetApplicationContext(ts, &honee)); 66 PetscCall(DMGetDimension(honee->dm, &dim)); 67 ceed = honee->ceed; 68 69 PetscCall(PetscNew(&monitor_ctx)); 70 PetscCall(HoneeCalculateDomainSize(honee, &monitor_ctx->volume)); 71 72 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(honee->elem_restr_q, &num_comp_q)); 73 PetscCall(QDataGet(ceed, honee->dm, domain_label, label_value, honee->elem_restr_x, honee->basis_x, honee->x_coord, &elem_restr_qd, &q_data, 74 &q_data_size)); 75 76 { // Get restriction and basis from the RHS function 77 CeedOperator *sub_ops, main_op = honee->op_ifunction ? honee->op_ifunction : honee->op_rhs_ctx->op; 78 CeedOperatorField op_field; 79 PetscInt sub_op_index = 0; // will be 0 for the volume op 80 81 PetscCallCeed(ceed, CeedOperatorCompositeGetSubList(main_op, &sub_ops)); 82 PetscCallCeed(ceed, CeedOperatorGetFieldByName(sub_ops[sub_op_index], "q", &op_field)); 83 84 PetscCallCeed(ceed, CeedOperatorFieldGetData(op_field, NULL, &elem_restr_q, &basis_q, NULL)); 85 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &newt_qfctx)); 86 } 87 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_totalke, basis_q, elem_restr_q, &elem_restr_totalke)); 88 89 switch (honee->phys->state_var) { 90 case STATEVAR_PRIMITIVE: 91 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Prim, MonitorTotalKineticEnergy_Prim_loc, &qf_monitor)); 92 break; 93 case STATEVAR_CONSERVATIVE: 94 PetscCallCeed(ceed, 95 CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Conserv, MonitorTotalKineticEnergy_Conserv_loc, &qf_monitor)); 96 break; 97 case STATEVAR_ENTROPY: 98 PetscCallCeed(ceed, 99 CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Entropy, MonitorTotalKineticEnergy_Entropy_loc, &qf_monitor)); 100 break; 101 } 102 103 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_monitor, newt_qfctx)); 104 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q", num_comp_q, CEED_EVAL_INTERP)); 105 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD)); 106 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q_data", q_data_size, CEED_EVAL_NONE)); 107 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_monitor, "v", num_comp_totalke, CEED_EVAL_NONE)); 108 109 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_monitor, NULL, NULL, &op_monitor)); 110 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 111 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE)); 112 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data)); 113 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "v", elem_restr_totalke, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 114 115 PetscCall(OperatorApplyContextCreate(honee->dm, NULL, ceed, op_monitor, NULL, NULL, NULL, NULL, &monitor_ctx->op_monitor_ctx)); 116 117 PetscCall(CeedOperatorCreateLocalVecs(op_monitor, DMReturnVecType(honee->dm), PETSC_COMM_SELF, NULL, &monitor_ctx->values)); 118 PetscCall(VecSetBlockSize(monitor_ctx->values, num_comp_totalke)); 119 monitor_ctx->num_comps = num_comp_totalke; 120 PetscCall(PetscMalloc1(num_comp_totalke, &monitor_ctx->sum_values)); 121 PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level)); 122 monitor_ctx->tab_level = tab_level + 1; 123 124 ctx->data = monitor_ctx; 125 ctx->data_destroy = (PetscCtxDestroyFn *)MonitorTotalKEDestroy; 126 127 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&newt_qfctx)); 128 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd)); 129 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q)); 130 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_totalke)); 131 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q)); 132 PetscCallCeed(ceed, CeedVectorDestroy(&q_data)); 133 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_monitor)); 134 PetscCallCeed(ceed, CeedOperatorDestroy(&op_monitor)); 135 PetscFunctionReturn(PETSC_SUCCESS); 136 } 137 138 PetscErrorCode TSMonitor_TotalKineticEnergy(TS ts, PetscInt step, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) { 139 MonitorTotalKE monitor_ctx = (MonitorTotalKE)ctx->data; 140 Honee honee; 141 MPI_Comm comm; 142 PetscMPIInt rank; 143 TSConvergedReason reason; 144 static const char *field_names[] = {"TotalKineticEnergy", "StrainDissipation", "DivergenceDissipation", "VolumeExpansion", "MuVorticity2"}; 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 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 152 153 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time)); 154 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc)); 155 156 PetscCall(ApplyCeedOperatorLocalToLocal(honee->Q_loc, monitor_ctx->values, monitor_ctx->op_monitor_ctx)); 157 158 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) { 159 PetscCall(VecStrideSum(monitor_ctx->values, i, &monitor_ctx->sum_values[i])); 160 } 161 162 if (rank == 0) PetscCallMPI(MPI_Reduce(MPI_IN_PLACE, monitor_ctx->sum_values, monitor_ctx->num_comps, MPIU_SCALAR, MPI_SUM, 0, comm)); 163 else PetscCallMPI(MPI_Reduce(monitor_ctx->sum_values, monitor_ctx->sum_values, monitor_ctx->num_comps, MPIU_SCALAR, MPI_SUM, 0, comm)); 164 165 if (rank == 0) 166 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) monitor_ctx->sum_values[i] /= monitor_ctx->volume; 167 168 if (ctx->format == PETSC_VIEWER_ASCII_CSV) { 169 if (!monitor_ctx->is_header_written) { 170 char buf[PETSC_MAX_PATH_LEN]; 171 const char *buf_const; 172 int major, minor, patch; 173 time_t t = time(NULL); 174 175 PetscCall(HoneeGetVersion(&major, &minor, &patch, NULL)); 176 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# honee_version: %d.%d.%d\n", major, minor, patch)); 177 PetscCall(HoneeGetGitVersion(&buf_const)); 178 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# honee_git_commit: %s\n", buf_const)); 179 PetscCall(PetscGetVersion(buf, sizeof(buf))); 180 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_version: %s\n", buf)); 181 PetscCall(PetscGetPetscDir(&buf_const)); 182 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_directory: %s\n", buf_const)); 183 PetscCall(PetscGetArchType(buf, sizeof(buf))); 184 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_arch: %s\n", buf)); 185 PetscCheck(strftime(buf, sizeof(buf), "%FT%T%z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed"); 186 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date: %s\n", buf)); 187 PetscCheck(strftime(buf, sizeof(buf), "%Z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed"); 188 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date_timezone: %s\n", buf)); 189 PetscCall(PetscGetUserName(buf, sizeof(buf))); 190 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# username: %s\n", buf)); 191 PetscCall(PetscGetHostName(buf, sizeof(buf))); 192 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# hostname: %s\n", buf)); 193 PetscCall(PetscGetWorkingDirectory(buf, sizeof(buf))); 194 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# working_directory: %s\n", buf)); 195 PetscCall(PetscViewerFileGetName(ctx->viewer, &buf_const)); 196 PetscCall(PetscGetFullPath(buf_const, buf, sizeof(buf))); 197 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# original_file_path: %s\n", buf)); 198 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "#\n")); 199 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "step,time,")); 200 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) { 201 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%s", field_names[i])); 202 if (i < monitor_ctx->num_comps - 1) PetscCall(PetscViewerASCIIPrintf(ctx->viewer, ",")); 203 } 204 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n")); 205 monitor_ctx->is_header_written = PETSC_TRUE; 206 } 207 208 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%" PetscInt_FMT ",%0.17e,", step, solution_time)); 209 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) { 210 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%0.17e", monitor_ctx->sum_values[i])); 211 if (i < monitor_ctx->num_comps - 1) PetscCall(PetscViewerASCIIPrintf(ctx->viewer, ",")); 212 } 213 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n")); 214 } else { 215 PetscCall(PetscViewerASCIIAddTab(ctx->viewer, monitor_ctx->tab_level)); 216 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) { 217 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%s: %0.12e ", field_names[i], monitor_ctx->sum_values[i])); 218 if (i == 0) PetscCall(PetscViewerASCIIUseTabs(ctx->viewer, PETSC_FALSE)); 219 } 220 PetscScalar sum_rates = 0; 221 for (PetscInt i = 1; i < monitor_ctx->num_comps - 1; i++) sum_rates += monitor_ctx->sum_values[i]; 222 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "SumOfRates: %0.12e ", sum_rates)); 223 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n")); 224 PetscCall(PetscViewerASCIIUseTabs(ctx->viewer, PETSC_TRUE)); 225 PetscCall(PetscViewerASCIISubtractTab(ctx->viewer, monitor_ctx->tab_level)); 226 } 227 PetscFunctionReturn(PETSC_SUCCESS); 228 } 229