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
MonitorCflCtxDestroy(MonitorCflCtx * monitor_ctx)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
SetupMontiorCfl(TS ts,PetscViewerAndFormat * ctx)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, &elem_restr_qd, &q_data, &q_data_size));
53 PetscCall(DMPlexCeedElemRestrictionQDataCreate(ceed, honee->dm, DMLABEL_DEFAULT, DMLABEL_DEFAULT_VALUE, 0, num_comp_cfl, &elem_restr_cfl));
54
55 { // Get restriction and basis from the RHS function
56 CeedOperator *sub_ops, main_op = honee->op_ifunction ? honee->op_ifunction : honee->op_rhs_ctx->op;
57 PetscInt sub_op_index = 0; // will be 0 for the volume op
58
59 PetscCallCeed(ceed, CeedOperatorCompositeGetSubList(main_op, &sub_ops));
60 PetscCallCeed(ceed, CeedOperatorGetContext(sub_ops[sub_op_index], &newt_qfctx));
61 }
62
63 switch (dim) {
64 case 2:
65 switch (honee->phys->state_var) {
66 case STATEVAR_PRIMITIVE:
67 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Prim, MonitorCFL_2D_Prim_loc, &qf_monitor));
68 break;
69 case STATEVAR_CONSERVATIVE:
70 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Conserv, MonitorCFL_2D_Conserv_loc, &qf_monitor));
71 break;
72 case STATEVAR_ENTROPY:
73 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_2D_Entropy, MonitorCFL_2D_Entropy_loc, &qf_monitor));
74 break;
75 }
76 break;
77 case 3:
78 switch (honee->phys->state_var) {
79 case STATEVAR_PRIMITIVE:
80 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Prim, MonitorCFL_3D_Prim_loc, &qf_monitor));
81 break;
82 case STATEVAR_CONSERVATIVE:
83 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Conserv, MonitorCFL_3D_Conserv_loc, &qf_monitor));
84 break;
85 case STATEVAR_ENTROPY:
86 PetscCallCeed(ceed, CeedQFunctionCreateInterior(ceed, 1, MonitorCFL_3D_Entropy, MonitorCFL_3D_Entropy_loc, &qf_monitor));
87 break;
88 }
89 break;
90 }
91 PetscCheck(qf_monitor, PetscObjectComm((PetscObject)ts), PETSC_ERR_SUP,
92 "Could not create CFL monitor QFunction for dim %" PetscInt_FMT " and state variable %s", dim, StateVariables[honee->phys->state_var]);
93
94 PetscBool is_advection;
95 PetscCall(PetscStrcmp("advection", honee->app_ctx->problem_name, &is_advection));
96 if (is_advection) {
97 NewtonianIdealGasContext newt_ctx;
98
99 PetscCall(PetscNew(&newt_ctx));
100 PetscCallCeed(ceed, CeedQFunctionContextDestroy(&newt_qfctx));
101 PetscCallCeed(ceed, CeedQFunctionContextCreate(ceed, &newt_qfctx));
102 PetscCallCeed(ceed, CeedQFunctionContextSetData(newt_qfctx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*newt_ctx), newt_ctx));
103 PetscCallCeed(ceed, CeedQFunctionContextSetDataDestroy(newt_qfctx, CEED_MEM_HOST, FreeContextPetsc));
104 }
105
106 PetscCallCeed(ceed, CeedQFunctionSetContext(qf_monitor, newt_qfctx));
107 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q", num_comp_q, CEED_EVAL_INTERP));
108 PetscCallCeed(ceed, CeedQFunctionAddInput(qf_monitor, "q_data", q_data_size, CEED_EVAL_NONE));
109 PetscCallCeed(ceed, CeedQFunctionAddOutput(qf_monitor, "v", num_comp_cfl, CEED_EVAL_NONE));
110
111 PetscCallCeed(ceed, CeedOperatorCreate(ceed, qf_monitor, NULL, NULL, &op_monitor));
112 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q", elem_restr_q, basis_q, CEED_VECTOR_ACTIVE));
113 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "q_data", elem_restr_qd, CEED_BASIS_NONE, q_data));
114 PetscCallCeed(ceed, CeedOperatorSetField(op_monitor, "v", elem_restr_cfl, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE));
115
116 PetscCall(OperatorApplyContextCreate(honee->dm, NULL, ceed, op_monitor, NULL, NULL, NULL, NULL, &monitor_ctx->op_monitor_ctx));
117
118 PetscCall(CeedOperatorCreateLocalVecs(op_monitor, DMReturnVecType(honee->dm), PETSC_COMM_SELF, NULL, &monitor_ctx->values));
119 PetscCall(VecSetBlockSize(monitor_ctx->values, num_comp_cfl));
120 monitor_ctx->num_comps = num_comp_cfl;
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 *)MonitorCflCtxDestroy;
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_cfl));
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
TSMonitor_Cfl(TS ts,PetscInt step,PetscReal solution_time,Vec Q,PetscViewerAndFormat * ctx)138 PetscErrorCode TSMonitor_Cfl(TS ts, PetscInt step, PetscReal solution_time, Vec Q, PetscViewerAndFormat *ctx) {
139 MonitorCflCtx monitor_ctx = (MonitorCflCtx)ctx->data;
140 Honee honee;
141 MPI_Comm comm;
142 TSConvergedReason reason;
143 PetscReal part_minmax[2], global_minmax[2], dt;
144
145 PetscFunctionBeginUser;
146 PetscCall(TSGetConvergedReason(ts, &reason));
147 if (!(step % ctx->view_interval == 0 || reason != TS_CONVERGED_ITERATING)) PetscFunctionReturn(PETSC_SUCCESS);
148 PetscCall(TSGetApplicationContext(ts, &honee));
149 PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
150
151 PetscCall(UpdateBoundaryValues(honee, honee->Q_loc, solution_time));
152 PetscCall(DMGlobalToLocal(honee->dm, Q, INSERT_VALUES, honee->Q_loc));
153
154 PetscCall(ApplyCeedOperatorLocalToLocal(honee->Q_loc, monitor_ctx->values, monitor_ctx->op_monitor_ctx));
155
156 PetscCall(VecMin(monitor_ctx->values, NULL, &part_minmax[0]));
157 PetscCall(VecMax(monitor_ctx->values, NULL, &part_minmax[1]));
158 // Need to get global min/max since `values` is only the local partition
159 PetscCall(PetscGlobalMinMaxReal(comm, part_minmax, global_minmax));
160 PetscCall(TSGetTimeStep(ts, &dt));
161 global_minmax[0] *= dt;
162 global_minmax[1] *= dt;
163
164 if (ctx->format == PETSC_VIEWER_ASCII_CSV) {
165 if (!monitor_ctx->is_header_written) {
166 char buf[PETSC_MAX_PATH_LEN];
167 const char *buf_const;
168 time_t t = time(NULL);
169
170 PetscCall(PetscGetVersion(buf, sizeof(buf)));
171 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_version: %s\n", buf));
172 PetscCall(PetscGetPetscDir(&buf_const));
173 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_directory: %s\n", buf_const));
174 PetscCall(PetscGetArchType(buf, sizeof(buf)));
175 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# petsc_arch: %s\n", buf));
176 PetscCheck(strftime(buf, sizeof(buf), "%FT%T%z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed");
177 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date: %s\n", buf));
178 PetscCheck(strftime(buf, sizeof(buf), "%Z", localtime(&t)), comm, PETSC_ERR_SYS, "strftime() call failed");
179 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# date_timezone: %s\n", buf));
180 PetscCall(PetscGetUserName(buf, sizeof(buf)));
181 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# username: %s\n", buf));
182 PetscCall(PetscGetHostName(buf, sizeof(buf)));
183 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# hostname: %s\n", buf));
184 PetscCall(PetscGetWorkingDirectory(buf, sizeof(buf)));
185 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# working_directory: %s\n", buf));
186
187 PetscCall(PetscViewerFileGetName(ctx->viewer, &buf_const));
188 PetscCall(PetscGetFullPath(buf_const, buf, sizeof(buf)));
189 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "# original_file_path: %s\n", buf));
190 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "#\n"));
191 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "step,time,mincfl,maxcfl\n"));
192 monitor_ctx->is_header_written = PETSC_TRUE;
193 }
194
195 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%" PetscInt_FMT ",%0.17e,", step, solution_time));
196 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "%0.17e,%0.17e", global_minmax[0], global_minmax[1]));
197 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "\n"));
198 } else {
199 PetscCall(PetscViewerASCIIAddTab(ctx->viewer, monitor_ctx->tab_level));
200 PetscCall(PetscViewerASCIIPrintf(ctx->viewer, "Min, Max CFL: %0.5g, %0.5g\n", global_minmax[0], global_minmax[1]));
201 PetscCall(PetscViewerASCIISubtractTab(ctx->viewer, monitor_ctx->tab_level));
202 }
203 PetscFunctionReturn(PETSC_SUCCESS);
204 }
205