From de1200e87205dff8d4e06c85f5aa6a631227c7aa Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Fri, 6 Dec 2024 11:15:32 +0100
Subject: [PATCH] Add adaptive mesh refinement and 2D test case
---
dart/src/wSolver.cpp | 75 +++++++++++++++++++++++++++
dart/src/wSolver.h | 1 +
dart/tests/amr.py | 117 +++++++++++++++++++++++++++++++++++++++++++
ext/amfe | 2 +-
4 files changed, 194 insertions(+), 1 deletion(-)
create mode 100644 dart/tests/amr.py
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index b9a7392..4e9f57f 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -206,6 +206,81 @@ void Solver::save(MshExport *mshWriter, std::string const &suffix)
std::cout << "done." << std::endl;
}
+/**
+ * @brief Compute refinment sensor and write new mesh sizes
+ * @todo Remesh internally
+ * @todo Investigate combination of gradient sensor with Laplacian sensor
+ * @todo Investigate effect of user-defined parameters (thresholds, refinement/coarsening factors, termination, etc.)
+ * @todo Add capability to define region of interest (LE, TE, exclude tip, etc.)
+ */
+bool Solver::refine(double minSize, tbox::MshExport *mshWriter, std::string const &suffix)
+{
+ // Compute deltaRho at cell in the streamwise direction
+ std::vector<double> dRhoE(pbl->msh->elems.size());
+ std::vector<double> dRho(pbl->msh->nodes.size()), mSize(pbl->msh->nodes.size());
+ tbb::parallel_for_each(pbl->fluid->adjMap.begin(), pbl->fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p) {
+ dRhoE[p.first->no - 1] = std::abs(pbl->fluid->rho->eval(*(p.first), phi, 0) - pbl->fluid->rho->eval(*(p.second[0]), phi, 0));
+ });
+ // Average deltaRho and compute mesh size at nodes
+ tbb::parallel_for_each(pbl->fluid->neMap.begin(), pbl->fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
+ dRho[p.first->row] = 0.;
+ mSize[p.first->row] = 0.;
+ for (auto e : p.second)
+ {
+ dRho[p.first->row] += dRhoE[e->no - 1];
+ if (pbl->nDim == 2)
+ mSize[p.first->row] += std::sqrt(2 / std::sqrt(3) * e->getDetJ(0));
+ else
+ mSize[p.first->row] += std::cbrt(2 / std::sqrt(2) * e->getDetJ(0));
+ }
+ dRho[p.first->row] /= p.second.size();
+ mSize[p.first->row] /= p.second.size();
+ });
+
+ // Sort the mesh size of the element according the their deltaRho
+ std::vector<std::pair<int, double>> sorter(pbl->msh->nodes.size());
+ tbb::parallel_for_each(pbl->msh->nodes.begin(), pbl->msh->nodes.end(), [&](Node *n) {
+ sorter[n->row] = std::pair<int, double>(n->row, dRho[n->row]);
+ });
+ std::sort(sorter.begin(), sorter.end(), [&](const std::pair<int, double> &a, const std::pair<int, double> &b) {
+ return a.second > b.second;
+ });
+ // Coarsen the 30% elements which have lowest dRho, and refine the 30% that have highest dRho
+ size_t refBnd = static_cast<size_t>(std::floor(0.3 * sorter.size()));
+ size_t crsBnd = static_cast<size_t>(std::floor(0.7 * sorter.size()));
+ tbb::parallel_for(tbb::blocked_range<size_t>(0, sorter.size()), [&](tbb::blocked_range<size_t> rng) {
+ for (size_t i = rng.begin(); i < rng.end(); ++i)
+ {
+ // Manually prevent refine in certain regions
+ // Node *n = pbl->msh->nodes[sorter[i].first];
+ // if (n->pos(1) > 0.98 * 14 && n->pos(1) < 1.02 * 14)
+ // if (mSize[sorter[i].first] < 0.045)
+ // continue;
+ if (i < refBnd)
+ mSize[sorter[i].first] *= 0.5;
+ else if (i > crsBnd)
+ mSize[sorter[i].first] *= 2.0;
+ else
+ mSize[sorter[i].first] *= 1.0;
+ }
+ });
+
+ // Write sensor and new mesh sizes
+ Results rSens;
+ rSens.scalars_at_nodes["sensor"] = &dRho;
+ mshWriter->saveList(pbl->msh->name + "_sensor" + suffix, rSens);
+ Results rSize;
+ rSize.scalars_at_nodes["size"] = &mSize;
+ mshWriter->saveList(pbl->msh->name + "_size" + suffix, rSize);
+
+ // Stop the refinement loop if the mesh size is lower that the specified target
+ std::cout << std::fixed << std::setprecision(5) << "min. size: " << *std::min_element(mSize.begin(), mSize.end()) << std::endl;
+ if (*std::min_element(mSize.begin(), mSize.end()) < minSize)
+ return true;
+ else
+ return false;
+}
+
/**
* @brief Compute total aerodynamic load coefficients on boundaries
*/
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index a988bca..ca26877 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -79,6 +79,7 @@ public:
void reset();
virtual Status run();
void save(tbox::MshExport *mshWriter, std::string const &suffix = "");
+ bool refine(double minSize, tbox::MshExport *mshWriter, std::string const &suffix = "");
int getRow(size_t i) const { return rows[i]; };
diff --git a/dart/tests/amr.py b/dart/tests/amr.py
new file mode 100644
index 0000000..baf7989
--- /dev/null
+++ b/dart/tests/amr.py
@@ -0,0 +1,117 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute transonic potential flow around a RAE 2822
+# Adrien Crovato
+#
+# Test the adaptive mesh refinement capability
+#
+# CAUTION
+# This test is not stable. The adaptive mesh refinement procedure must be improved.
+
+import dart.utils as floU
+import dart.default as floD
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+import math
+
+def main():
+ # timer
+ tms = fwk.Timers()
+ tms['total'].start()
+
+ # define flow variables
+ alpha = 2*math.pi/180 # must be zero for this testfile
+ M_inf = 0.72
+ c_ref = 1
+ dim = 2
+
+ # define initial mesh parameters
+ mName = 'rae2822'
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 2.0, 'msTe' : 0.1, 'msLe' : 0.1}
+
+ # AMR loop
+ counter = 0
+ Cls = []
+ Cds = []
+ Cms = []
+ while True:
+ # mesh the geometry
+ tms['msh'].start()
+ msh, gmshWriter = floD.mesh(dim, f'models/{mName}.geo', pars, ['field', 'airfoil', 'downstream'])
+ tms['msh'].stop()
+
+ # set the problem and add fluid, initial/boundary conditions
+ tms['pre'].start()
+ pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil')
+ tms['pre'].stop()
+
+ # solve problem
+ tms['solver'].start()
+ solver = floD.newton(pbl)
+ solver.run()
+ solver.save(gmshWriter)
+ stop = solver.refine(0.005, gmshWriter)
+ tms['solver'].stop()
+
+ # extract data
+ Cls.append(solver.Cl)
+ Cds.append(solver.Cd)
+ Cms.append(solver.Cm)
+ Cp = floU.extract(bnd.groups[0].tag.elems, solver.Cp)
+ tboxU.write(Cp, f'Cp_airfoil_{counter}.dat', '%1.4e', ',', 'x, y, cp', '')
+
+ # stop criterion
+ if counter > 10 or stop:
+ break
+ else:
+ counter += 1
+ pars['-bgm'] = f'{mName}_size.pos'
+ pars['counter'] = counter # pass the counter to force remesh
+
+ # display results
+ print('{0:>8s} {1:>8s} {2:>8s} {3:>8s} {4:>8s} {5:>8s}'.format('run', 'Mach', 'AoA', 'Cl', 'Cd', 'Cm'))
+ i = 0
+ while i < len(Cls):
+ print('{0:8d} {1:8.3f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f}'.format(i, M_inf, alpha*180/math.pi, Cls[i], Cds[i], Cms[i]))
+ i = i+1
+ print('\n')
+
+ # display timers
+ tms['total'].stop()
+ print('CPU statistics')
+ print(tms)
+
+ # visualize solution and plot results
+ floD.initViewer(pbl)
+ floD.plot(Cp[:,0], Cp[:,-1], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+
+ # check results
+ tests = CTests()
+ tests.add(CTest('iteration count', solver.nIt, 15, 3, forceabs=True))
+ tests.add(CTest('refinement count', counter, 4, 1, forceabs=True))
+ tests.add(CTest('Cl', solver.Cl, 0.85, 1e-2))
+ tests.add(CTest('Cd', solver.Cd, 0.0025, 1e-2))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == "__main__":
+ main()
diff --git a/ext/amfe b/ext/amfe
index c11ff6e..fb9deb4 160000
--- a/ext/amfe
+++ b/ext/amfe
@@ -1 +1 @@
-Subproject commit c11ff6ee646765a19f9d43496d354e939a84b408
+Subproject commit fb9deb4bef7e323b709568d214205a15dd78359d
--
GitLab