Skip to content

Commit b40e595

Browse files
Bot-Enigma-0Bot-Enigma-0
committed
apply suggestions
1 parent 39a2375 commit b40e595

File tree

3 files changed

+205
-117
lines changed

3 files changed

+205
-117
lines changed

‎SU2_CFD/include/drivers/CSinglezoneDriver.hpp‎

Lines changed: 9 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -46,52 +46,19 @@ class CSinglezoneDriver : public CDriver {
4646
* \return Boolean indicating whether the problem is converged.
4747
*/
4848
virtual bool GetTimeConvergence() const;
49-
vector<passivedouble> v_estimate;
49+
su2matrix<passivedouble> v_estimate;
5050

5151
/*!
52-
* \brief Seed derivatives for all solvers in the zone.
53-
* \param[in] derivatives - Vector of derivative values
54-
*/
55-
void SeedAllDerivatives(const vector<passivedouble>& derivatives, CVariable* nodes, CGeometry *geometry) {
56-
unsigned long offset = 0;
57-
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
58-
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
59-
if (!solver) continue;
60-
61-
unsigned short nVar = solver->GetnVar();
62-
unsigned long nPoint = geometry->GetnPointDomain();
63-
unsigned long nTotal = nVar * nPoint;
64-
65-
vector<passivedouble> solver_derivs(derivatives.begin() + offset,
66-
derivatives.begin() + offset + nTotal);
67-
solver->SetSolutionDerivatives(solver_derivs, nodes);
68-
offset += nTotal;
69-
}
70-
}
52+
* \brief Seed derivatives for all solvers in the zone.
53+
* \param[in] derivatives - Matrix of derivative values
54+
*/
55+
void SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry *geometry);
7156

7257
/*!
73-
* \brief Extract derivatives from all solvers in the zone.
74-
* \param[out] derivatives - Vector to store derivative values
75-
*/
76-
void GetAllDerivatives(vector<passivedouble>& derivatives, CVariable* nodes, CGeometry *geometry) {
77-
unsigned long offset = 0;
78-
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
79-
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
80-
if (!solver) continue;
81-
82-
unsigned short nVar = solver->GetnVar();
83-
unsigned long nPoint = geometry->GetnPointDomain();
84-
unsigned long nTotal = nVar * nPoint;
85-
86-
vector<passivedouble> solver_derivs(nTotal);
87-
solver->GetSolutionDerivatives(solver_derivs, nodes);
88-
89-
for (unsigned long i = 0; i < nTotal; i++) {
90-
derivatives[offset + i] = solver_derivs[i];
91-
}
92-
offset += nTotal;
93-
}
94-
}
58+
* \brief Extract derivatives from all solvers in the zone.
59+
* \param[out] derivatives - Matrix to store derivative values
60+
*/
61+
void GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry *geometry);
9562

9663
public:
9764

‎SU2_CFD/include/solvers/CSolver.hpp‎

Lines changed: 0 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -4354,34 +4354,6 @@ inline void CustomSourceResidual(CGeometry *geometry, CSolver **solver_container
43544354
AD::EndNoSharedReading();
43554355
}
43564356

4357-
/*!
4358-
* \brief Seed derivatives for all solution variables using a flat vector.
4359-
* \param[in] derivatives - Flat vector of derivative values (size nVar * nPoint)
4360-
*/
4361-
void SetSolutionDerivatives(const vector<passivedouble>& derivatives, CVariable* nodes) {
4362-
unsigned long offset = 0;
4363-
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
4364-
su2double* solution = nodes->GetSolution(iPoint);
4365-
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
4366-
SU2_TYPE::SetDerivative(solution[iVar], derivatives[offset++]);
4367-
}
4368-
}
4369-
}
4370-
4371-
/*!
4372-
* \brief Get derivatives from all solution variables into a flat vector.
4373-
* \param[out] derivatives - Flat vector to store derivative values (size nVar * nPoint)
4374-
*/
4375-
void GetSolutionDerivatives(vector<passivedouble>& derivatives, CVariable* nodes) const {
4376-
unsigned long offset = 0;
4377-
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
4378-
su2double* solution = nodes->GetSolution(iPoint);
4379-
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
4380-
derivatives[offset++] = SU2_TYPE::GetDerivative(solution[iVar]);
4381-
}
4382-
}
4383-
}
4384-
43854357
protected:
43864358
/*!
43874359
* \brief Allocate the memory for the verification solution, if necessary.

‎SU2_CFD/src/drivers/CSinglezoneDriver.cpp‎

Lines changed: 196 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -313,82 +313,231 @@ bool CSinglezoneDriver::GetTimeConvergence() const{
313313
return output_container[ZONE_0]->GetCauchyCorrectedTimeConvergence(config_container[ZONE_0]);
314314
}
315315

316+
void CSinglezoneDriver::SeedAllDerivatives(const su2matrix<passivedouble>& derivatives, CGeometry *geometry) {
317+
unsigned long pointOffset = 0;
318+
unsigned long varOffset = 0;
319+
320+
if (rank == MASTER_NODE) {
321+
std::cout << "runnin SeedAllDerivatives function" << endl;
322+
std::cout << "matrix size: " << derivatives.rows() << " x " << derivatives.cols() << endl;
323+
}
324+
325+
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
326+
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
327+
if (!solver) continue;
328+
329+
CVariable* nodes = solver->GetNodes();
330+
if (!nodes) continue;
331+
332+
unsigned short nVar = solver->GetnVar();
333+
unsigned long nPoint = geometry->GetnPointDomain();
334+
335+
if (rank == MASTER_NODE) {
336+
std::cout << "solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << endl;
337+
}
338+
339+
// seed derivatives (removed function from csolver to make debug easier)
340+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
341+
su2double* solution = nodes->GetSolution(iPoint);
342+
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
343+
double derivative_value = derivatives(pointOffset + iPoint, varOffset + iVar);
344+
345+
// debugging value and derivative
346+
double original_value = SU2_TYPE::GetValue(solution[iVar]);
347+
SU2_TYPE::SetDerivative(solution[iVar], derivative_value);
348+
double extracted_derivative = SU2_TYPE::GetDerivative(solution[iVar]);
349+
350+
// see if forward mode works- output for the first few points -
351+
// WORKS WHEN RUNNING SU2_CFD_DIRECTDIFF, when SU2_CFD is run gives "Got= 0"
352+
if (iPoint < 5 && iVar == 0 && rank == MASTER_NODE) {
353+
std::cout << "Point " << iPoint << ", Var " << iVar << ": Set=" << derivative_value
354+
<< ", Got=" << extracted_derivative << ", Value=" << original_value << std::endl;
355+
}
356+
}
357+
}
358+
pointOffset += nPoint;
359+
varOffset += nVar;
360+
}
361+
362+
if (rank == MASTER_NODE) {
363+
std::cout << "Finished SeedAllDerivatives" << endl;
364+
}
365+
}
366+
367+
void CSinglezoneDriver::GetAllDerivatives(su2matrix<passivedouble>& derivatives, CGeometry *geometry) {
368+
unsigned long pointOffset = 0;
369+
unsigned long varOffset = 0;
370+
371+
if (rank == MASTER_NODE) {
372+
std::cout << "runnin GetAllDerivatives function" << endl;
373+
}
374+
375+
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
376+
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
377+
if (!solver) continue;
378+
379+
CVariable* nodes = solver->GetNodes();
380+
if (!nodes) continue;
381+
382+
unsigned short nVar = solver->GetnVar();
383+
unsigned long nPoint = geometry->GetnPointDomain();
384+
385+
if (rank == MASTER_NODE) {
386+
std::cout << "solver " << iSol << " has nVar=" << nVar << ", nPoint=" << nPoint << endl;
387+
}
388+
389+
// extract derivatives for solver
390+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
391+
su2double* solution = nodes->GetSolution(iPoint);
392+
for (unsigned short iVar = 0; iVar < nVar; iVar++) {
393+
double derivative_value = SU2_TYPE::GetDerivative(solution[iVar]);
394+
derivatives(pointOffset + iPoint, varOffset + iVar) = derivative_value;
395+
396+
// Debug: check if derivatives are being extracted correctly
397+
if (iPoint == 0 && iVar == 0 && rank == MASTER_NODE) {
398+
std::cout << "extracted derivative for solver " << iSol << ", point 0, var 0: " << derivative_value << std::endl;
399+
}
400+
}
401+
}
402+
pointOffset += nPoint;
403+
varOffset += nVar;
404+
}
405+
406+
if (rank == MASTER_NODE) {
407+
std::cout << "Finished GetAllDerivatives" << endl;
408+
}
409+
}
410+
316411
void CSinglezoneDriver::PreRunSpectralRadius() {
412+
if (!config_container[ZONE_0]->GetSpectralRadius_Analysis()) return;
413+
317414
CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0];
318-
CVariable* nodes = nullptr;
319415

320-
// Get total number of variables
416+
// get total number of variables and points
321417
unsigned long nVarTotal = 0;
322418
unsigned long nPoint = geometry->GetnPointDomain();
323419
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
324-
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
325-
if (solver) nVarTotal += solver->GetnVar();
420+
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
421+
if (solver) nVarTotal += solver->GetnVar();
326422
}
327423

328-
unsigned long nTotal = nVarTotal * nPoint;
329-
330-
// Initialize power iteration vector
331-
if (v_estimate.empty()) {
332-
v_estimate.resize(nTotal);
333-
passivedouble norm = 0.0;
334-
for (unsigned long i = 0; i < nTotal; i++) {
335-
v_estimate[i] = static_cast<passivedouble>(rand()) / RAND_MAX;
336-
norm += v_estimate[i] * v_estimate[i];
424+
// initialize power iteration matrix if needed
425+
if (v_estimate.rows() != nPoint || v_estimate.cols() != nVarTotal) {
426+
v_estimate.resize(nPoint, nVarTotal);
427+
428+
passivedouble norm = 0.0;
429+
// initialize with random values
430+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
431+
for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
432+
v_estimate(iPoint, iVar) = static_cast<passivedouble>(rand()) / RAND_MAX;
433+
norm += v_estimate(iPoint, iVar) * v_estimate(iPoint, iVar);
434+
}
435+
}
436+
437+
// norm calculation
438+
passivedouble global_norm;
439+
SU2_MPI::Allreduce(&norm, &global_norm, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
440+
global_norm = sqrt(global_norm);
441+
442+
// normalize the matrix
443+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
444+
for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
445+
v_estimate(iPoint, iVar) /= global_norm;
337446
}
338-
norm = sqrt(norm);
339-
for (unsigned long i = 0; i < nTotal; i++) v_estimate[i] /= norm;
447+
}
340448
}
341449

342-
// Seed derivatives for all solvers
343-
SeedAllDerivatives(v_estimate, nodes, geometry);
344-
if (rank == MASTER_NODE) {std::cout << "seeding done" << endl;}; //see where code fails
450+
// seed all solvers
451+
SeedAllDerivatives(v_estimate, geometry);
452+
if (rank == MASTER_NODE) { std::cout << "seeding done" << endl; }
345453
}
346454

347455
void CSinglezoneDriver::PostRunSpectralRadius() {
348-
if (rank == MASTER_NODE) {std::cout << "running post spectral function" << endl;}; //see where code fails
456+
if (!config_container[ZONE_0]->GetSpectralRadius_Analysis()) return;
457+
458+
if (rank == MASTER_NODE) {
459+
std::cout << "Running PostRunSpectralRadius" << endl;
460+
}
349461

350462
CGeometry* geometry = geometry_container[ZONE_0][INST_0][MESH_0];
351-
CVariable* nodes = nullptr;
352463

353-
// Get total number of variables
464+
// get total number of variables and points
354465
unsigned long nVarTotal = 0;
355466
unsigned long nPoint = geometry->GetnPointDomain();
356467
for (auto iSol = 0u; iSol < MAX_SOLS; ++iSol) {
357-
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
358-
if (solver) nVarTotal += solver->GetnVar();
468+
CSolver* solver = solver_container[ZONE_0][INST_0][MESH_0][iSol];
469+
if (solver) nVarTotal += solver->GetnVar();
470+
}
471+
472+
if (rank == MASTER_NODE) {
473+
std::cout << "total var: " << nVarTotal << ", total points: " << nPoint << endl;
359474
}
360-
unsigned long nTotal = nVarTotal * nPoint;
361-
if (rank == MASTER_NODE) {std::cout <<"extracting tangents" << endl;}; //see where code fails
362475

363-
// Extract derivatives from all solvers
364-
vector<passivedouble> w(nTotal);
365-
GetAllDerivatives(w, nodes, geometry);
366-
if (rank == MASTER_NODE) {std::cout << "computing eigen" << endl;}; //see where code fails
476+
// extract derivatives from all solvers
477+
su2matrix<passivedouble> w(nPoint, nVarTotal);
478+
479+
if (rank == MASTER_NODE) {
480+
std::cout << "matrix w initialized with size: " << w.rows() << " x " << w.cols() << endl;
481+
}
482+
483+
GetAllDerivatives(w, geometry);
484+
485+
if (rank == MASTER_NODE) {
486+
std::cout << "derivatives extracted" << endl;
487+
}
367488

368-
// Compute norm and update eigen estimate
489+
// compute norm
369490
passivedouble rho_new = 0.0;
370-
for (unsigned long i = 0; i < nTotal; i++) rho_new += w[i] * w[i];
371-
rho_new = sqrt(rho_new);
491+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
492+
for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
493+
rho_new += w(iPoint, iVar) * w(iPoint, iVar);
494+
}
495+
}
496+
497+
passivedouble global_rho_new;
498+
SU2_MPI::Allreduce(&rho_new, &global_rho_new, 1, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
499+
global_rho_new = sqrt(global_rho_new);
500+
501+
if (rank == MASTER_NODE) {
502+
std::cout << "Global rho_new: " << global_rho_new << endl;
503+
}
372504

373-
for (unsigned long i = 0; i < nTotal; i++) v_estimate[i] = w[i] / rho_new;
505+
// update eigenvector estimate
506+
for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
507+
for (unsigned short iVar = 0; iVar < nVarTotal; iVar++) {
508+
v_estimate(iPoint, iVar) = w(iPoint, iVar) / global_rho_new;
509+
}
510+
}
374511

375-
// Check convergence
512+
// check convergence
376513
static passivedouble rho_old = 0.0;
377514
static int iter = 0;
378515

379-
int max_iter = 500;
380-
passivedouble tol = 1e-7;
516+
int max_iter = config_container[ZONE_0]->GetnPowerMethodIter();
517+
su2double tol = config_container[ZONE_0]->GetPowerMethodTolerance();
381518

382-
if (abs(rho_new - rho_old) < tol || iter >= max_iter) {
383-
if (rank == MASTER_NODE) {
384-
std::cout << "Spectral Radius Estimate: " << rho_new << " after " << iter << " iterations." << std::endl;
385-
}
386-
// Reset for next analysis if needed
387-
rho_old = 0.0;
388-
iter = 0;
389-
v_estimate.clear();
519+
if (rank == MASTER_NODE) {
520+
std::cout << "Power Iterations requested = " << max_iter << endl;
521+
std::cout << "Power tolerance requested = " << tol << endl;
522+
std::cout << "Current iteration = " << iter << endl;
523+
std::cout << "Current rho_old = " << rho_old << endl;
524+
std::cout << "Current global_rho_new = " << global_rho_new << endl;
525+
}
526+
527+
if (abs(global_rho_new - rho_old) < tol || iter >= max_iter) {
528+
if (rank == MASTER_NODE) {
529+
std::cout << "Spectral Radius Estimate: " << global_rho_new << " after " << iter << " iterations." << std::endl;
530+
}
531+
// Reset
532+
rho_old = 0.0;
533+
iter = 0;
534+
v_estimate.resize(0, 0);
390535
} else {
391-
rho_old = rho_new;
392-
iter++;
536+
rho_old = global_rho_new;
537+
iter++;
538+
539+
if (rank == MASTER_NODE) {
540+
std::cout << "Continuing to iteration " << iter << " with rho_old = " << rho_old << endl;
541+
}
393542
}
394-
}
543+
}

0 commit comments

Comments
 (0)