Add adaptive mesh refinement and 2D test case

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
  */

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]; };
 

dart/tests/amr.py

+#!/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()

amfe @ fb9deb4b

-Subproject commit c11ff6ee646765a19f9d43496d354e939a84b408
+Subproject commit fb9deb4bef7e323b709568d214205a15dd78359d

[Adrien Crovato]

The AMR feature works quite well in 2D, but not in 3D. Several aspects must be considered and improved:

• The remeshing process should be handled internally and not through Gsmh in order to decrease the computational time.
• The mesh should be adapted by splitting or agglomerating the cells rather than being entirely regenerated from scratch.
• The refinement sensor is currently based on the difference between the density of the current and the upwind cell. A sensor based on the density gradient and/or on the Laplacian should be considered. Alternatively, the (adjoint) mesh sensitives could be used as a sensor.
• The refinement process should be limited to certain regions of interest.
• The refinement/coarsening process relies on user-defined parameters (refinement/coarsening thresholds and factors, termination criterion, etc.) which depend on the case being solved. This should be investigated to find an optimal set of parameters for standard cases.
• A restart capability should be implemented.

As a final note, it is unlikely that AMR will be faster than a single calculation performed directly on a "good" grid. However, it could alleviate the burden of generating this "good" grid from scratch and give insights into the flow pattern.