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
MonitorTotalKEDestroy(MonitorTotalKE * monitor_ctx)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 */
CreateElemRestrColloc_CompMajor(Ceed ceed,CeedInt num_comp,CeedBasis basis,CeedElemRestriction elem_restr_base,CeedElemRestriction * elem_restr_collocated)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
SetupMontiorTotalKineticEnergy(TS ts,PetscViewerAndFormat * ctx)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 PetscInt num_comp_totalke = 5, dim, tab_level;
57 MonitorTotalKE monitor_ctx;
58 CeedQFunctionContext newt_qfctx;
59 PetscBool is_ascii;
60
61 PetscFunctionBeginUser;
62 PetscCall(PetscObjectTypeCompare((PetscObject)ctx->viewer, PETSCVIEWERASCII, &is_ascii));
63 PetscCheck(is_ascii, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP, "Only supports ASCII viewers");
64 PetscCall(TSGetApplicationContext(ts, &honee));
65 PetscCall(DMGetDimension(honee->dm, &dim));
66 ceed = honee->ceed;
67
68 PetscCall(PetscNew(&monitor_ctx));
69 PetscCall(HoneeCalculateDomainSize(honee, &monitor_ctx->volume));
70 { // Get restriction and basis from the RHS function
71 CeedOperator *sub_ops, main_op = honee->op_ifunction ? honee->op_ifunction : honee->op_rhs_ctx->op;
72 PetscInt sub_op_index = 0; // will be 0 for the volume op
73
74 PetscCallCeed(ceed, CeedOperatorCompositeGetSubList(main_op, &sub_ops));
75 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &newt_qfctx));
76 }
77 PetscCall(QDataGet(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, &elem_restr_qd, &q_data, &q_data_size));
78 PetscCall(DMPlexCeedElemRestrictionCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &elem_restr_q));
79 PetscCall(DMPlexCeedBasisCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, 0, &basis_q));
80 PetscCall(CreateElemRestrColloc_CompMajor(ceed, num_comp_totalke, basis_q, elem_restr_q, &elem_restr_totalke));
81 PetscCallCeed(ceed, CeedElemRestrictionGetNumComponents(elem_restr_q, &num_comp_q));
82
83 switch (honee->phys->state_var) {
84 case STATEVAR_PRIMITIVE:
85 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Prim, MonitorTotalKineticEnergy_Prim_loc, &qf_monitor));
86 break;
87 case STATEVAR_CONSERVATIVE:
88 PetscCallCeed(ceed,
89 CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Conserv, MonitorTotalKineticEnergy_Conserv_loc, &qf_monitor));
90 break;
91 case STATEVAR_ENTROPY:
92 PetscCallCeed(ceed,
93 CeedQFunctionCreateInterior(ceed, 1, MonitorTotalKineticEnergy_Entropy, MonitorTotalKineticEnergy_Entropy_loc, &qf_monitor));
94 break;
95 }
96
97 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_monitor, newt_qfctx));
98 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q", num_comp_q, CEED_EVAL_INTERP));
99 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "Grad_q", num_comp_q * dim, CEED_EVAL_GRAD));
100 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q_data", q_data_size, CEED_EVAL_NONE));
101 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_monitor, "v", num_comp_totalke, CEED_EVAL_NONE));
102
103 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_monitor, NULL, NULL, &op_monitor));
104 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
105 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "Grad_q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
106 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data));
107 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "v", elem_restr_totalke, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
108
109 PetscCall(OperatorApplyContextCreate(honee->dm, NULL, ceed, op_monitor, NULL, NULL, NULL, NULL, &monitor_ctx->op_monitor_ctx));
110
111 PetscCall(CeedOperatorCreateLocalVecs(op_monitor, DMReturnVecType(honee->dm), PETSC_COMM_SELF, NULL, &monitor_ctx->values));
112 PetscCall(VecSetBlockSize(monitor_ctx->values, num_comp_totalke));
113 monitor_ctx->num_comps = num_comp_totalke;
114 PetscCall(PetscMalloc1(num_comp_totalke, &monitor_ctx->sum_values));
115 PetscCall(PetscObjectGetTabLevel((PetscObject)ts, &tab_level));
116 monitor_ctx->tab_level = tab_level + 1;
117
118 ctx->data = monitor_ctx;
119 ctx->data_destroy = (PetscCtxDestroyFn *)MonitorTotalKEDestroy;
120
121 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&newt_qfctx));
122 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_qd));
123 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_q));
124 PetscCallCeed(ceed, CeedElemRestrictionDestroy(&elem_restr_totalke));
125 PetscCallCeed(ceed, CeedBasisDestroy(&basis_q));
126 PetscCallCeed(ceed, CeedVectorDestroy(&q_data));
127 PetscCallCeed(ceed, CeedQFunctionDestroy(&qf_monitor));
128 PetscCallCeed(ceed, CeedOperatorDestroy(&op_monitor));
129 PetscFunctionReturn(PETSC_SUCCESS);
130 }
131
TSMonitor_TotalKineticEnergy(TS ts,PetscInt step,PetscReal solution_time,Vec Q,PetscViewerAndFormat * ctx)132 PetscErrorCode TSMonitor_TotalKineticEnergy(TS ts, PetscInt step, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) {
133 MonitorTotalKE monitor_ctx = (MonitorTotalKE)ctx->data;
134 Honee honee;
135 MPI_Comm comm;
136 PetscMPIInt rank;
137 TSConvergedReason reason;
138 static const char *field_names[] = {"TotalKineticEnergy", "StrainDissipation", "DivergenceDissipation", "VolumeExpansion", "MuVorticity2"};
139
140 PetscFunctionBeginUser;
141 PetscCall(TSGetConvergedReason(ts, &reason));
142 if (!(step % ctx->view_interval == 0 || reason != TS_CONVERGED_ITERATING)) PetscFunctionReturn(PETSC_SUCCESS);
143 PetscCall(TSGetApplicationContext(ts, &honee));
144 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
145 PetscCallMPI(MPI_Comm_rank(comm, &rank));
146
147 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time));
148 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc));
149
150 PetscCall(ApplyCeedOperatorLocalToLocal(honee->Q_loc, monitor_ctx->values, monitor_ctx->op_monitor_ctx));
151
152 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) {
153 PetscCall(VecStrideSum(monitor_ctx->values, i, &monitor_ctx->sum_values[i]));
154 }
155
156 if (rank == 0) PetscCallMPI(MPI_Reduce(MPI_IN_PLACE, monitor_ctx->sum_values, monitor_ctx->num_comps, MPIU_SCALAR, MPI_SUM, 0, comm));
157 else PetscCallMPI(MPI_Reduce(monitor_ctx->sum_values, monitor_ctx->sum_values, monitor_ctx->num_comps, MPIU_SCALAR, MPI_SUM, 0, comm));
158
159 if (rank == 0)
160 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) monitor_ctx->sum_values[i] /= monitor_ctx->volume;
161
162 if (ctx->format == PETSC_VIEWER_ASCII_CSV) {
163 if (!monitor_ctx->is_header_written) {
164 char buf[PETSC_MAX_PATH_LEN];
165 const char *buf_const;
166 int major, minor, patch;
167 time_t t = time(NULL);
168
169 PetscCall(HoneeGetVersion(&major, &minor, &patch, NULL));
170 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# honee_version: %d.%d.%d\n", major, minor, patch));
171 PetscCall(HoneeGetGitVersion(&buf_const));
172 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# honee_git_commit: %s\n", buf_const));
173 PetscCall(PetscGetVersion(buf, sizeof(buf)));
174 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_version: %s\n", buf));
175 PetscCall(PetscGetPetscDir(&buf_const));
176 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_directory: %s\n", buf_const));
177 PetscCall(PetscGetArchType(buf, sizeof(buf)));
178 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_arch: %s\n", buf));
179 PetscCheck(strftime(buf, sizeof(buf), "%FT%T%z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed");
180 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date: %s\n", buf));
181 PetscCheck(strftime(buf, sizeof(buf), "%Z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed");
182 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date_timezone: %s\n", buf));
183 PetscCall(PetscGetUserName(buf, sizeof(buf)));
184 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# username: %s\n", buf));
185 PetscCall(PetscGetHostName(buf, sizeof(buf)));
186 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# hostname: %s\n", buf));
187 PetscCall(PetscGetWorkingDirectory(buf, sizeof(buf)));
188 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# working_directory: %s\n", buf));
189 PetscCall(PetscViewerFileGetName(ctx->viewer, &buf_const));
190 PetscCall(PetscGetFullPath(buf_const, buf, sizeof(buf)));
191 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# original_file_path: %s\n", buf));
192 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "#\n"));
193 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "step,time,"));
194 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) {
195 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%s", field_names[i]));
196 if (i < monitor_ctx->num_comps - 1) PetscCall(PetscViewerASCIIPrintf(ctx->viewer, ","));
197 }
198 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n"));
199 monitor_ctx->is_header_written = PETSC_TRUE;
200 }
201
202 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%" PetscInt_FMT ",%0.17e,", step, solution_time));
203 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) {
204 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%0.17e", monitor_ctx->sum_values[i]));
205 if (i < monitor_ctx->num_comps - 1) PetscCall(PetscViewerASCIIPrintf(ctx->viewer, ","));
206 }
207 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n"));
208 } else {
209 PetscCall(PetscViewerASCIIAddTab(ctx->viewer, monitor_ctx->tab_level));
210 for (PetscInt i = 0; i < monitor_ctx->num_comps; i++) {
211 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%s: %0.12e ", field_names[i], monitor_ctx->sum_values[i]));
212 if (i == 0) PetscCall(PetscViewerASCIIUseTabs(ctx->viewer, PETSC_FALSE));
213 }
214 PetscScalar sum_rates = 0;
215 for (PetscInt i = 1; i < monitor_ctx->num_comps - 1; i++) sum_rates += monitor_ctx->sum_values[i];
216 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "SumOfRates: %0.12e ", sum_rates));
217 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n"));
218 PetscCall(PetscViewerASCIIUseTabs(ctx->viewer, PETSC_TRUE));
219 PetscCall(PetscViewerASCIISubtractTab(ctx->viewer, monitor_ctx->tab_level));
220 }
221 PetscFunctionReturn(PETSC_SUCCESS);
222 }
223