feat(project): main init commit

build-config/CLI11/meson.build

@@ -0,0 +1,2 @@
+cli11_proj = subproject('cli11')
+cli11_dep = cli11_proj.get_variable('CLI11_dep')

build-config/libcomposition/meson.build

@@ -0,0 +1,13 @@
+composition_p = subproject('libcomposition',
+  default_options: [
+    'pkg_config=false',
+    'build_tests=false',
+    'build_examples=false'
+])
+comp_dep  = composition_p.get_variable('composition_dep')
+libcomposition = composition_p.get_variable('libcomposition')
+spw_dep = composition_p.get_variable('species_weight_dep')
+composition_dep = [
+    comp_dep,
+    spw_dep,
+]

build-config/libconfig/meson.build

@@ -0,0 +1,8 @@
+config_p = subproject('libconfig',
+default_options:[
+  'default_library=static',
+  'pkg_config=false',
+  'build_tests=false',
+  'build_examples=false'
+])
+config_dep  = config_p.get_variable('config_dep')

build-config/meson.build

@@ -0,0 +1,4 @@
+subdir('mfem')
+subdir('libconfig')
+subdir('CLI11')
+subdir('libcomposition')

build-config/mfem/meson.build

@@ -0,0 +1,16 @@
+cmake = import('cmake')
+mfem_cmake_options = cmake.subproject_options()
+mfem_cmake_options.add_cmake_defines({
+    'MFEM_ENABLE_EXAMPLES': 'OFF',
+    'MFEM_ENABLE_TESTING': 'OFF',
+    'MFEM_ENABLE_MINIAPPS': 'OFF',
+    'MFEM_USE_BENCMARK': 'OFF',
+    'BUILD_SHARED_LIBS': 'OFF',
+    'BUILD_STATIC_LIBS': 'ON',
+})
+mfem_cmake_options.set_install(true)
+
+mfem_sp = cmake.subproject(
+    'mfem',
+    options: mfem_cmake_options)
+mfem_dep = mfem_sp.dependency('mfem')

default.toml

@@ -0,0 +1,10 @@
+[main]
+core_steepness = 1.0
+flattening = 0.2
+include_external_domain = false
+order = 3
+r_core = 1.5
+r_infinity = 6.0
+r_instability = 1e-14
+r_star = 5.0
+refinement_levels = 3

main.cpp

@@ -0,0 +1,187 @@
+#include "mfem.hpp"
+#include "fourdst/config/config.h"
+#include "fourdst/composition/composition.h"
+#include "CLI/CLI.hpp"
+
+#include <functional>
+
+struct Options {
+    std::string input_mesh = "stroid.mesh";
+    std::string vishost = "localhost";
+    int visport = 19916;
+};
+
+void ViewMesh(
+    const std::string& vishost,
+    int visport,
+    const mfem::Mesh& mesh,
+    mfem::GridFunction& gf,
+    const std::string& title,
+    bool first_step = true
+) {
+    mfem::socketstream sol_sock(vishost.c_str(), visport);
+    if (!sol_sock.is_open()) {
+        std::cerr << "Unable to connect to GLVis server at "
+                  << vishost << ':' << visport << std::endl;
+        std::exit(1);
+    }
+
+    sol_sock.precision(8);
+    sol_sock << "solution\n" << mesh << gf;
+    if (first_step) {
+        sol_sock << "window_title '" << title << "'\n";
+        sol_sock << "keys 'iIzzMaagpmtppc'";
+    }
+    sol_sock << std::flush;
+}
+
+double density(const mfem::Vector& x) {
+    double radius = x.Norml2();
+    constexpr double R = 5;
+    constexpr double rho0 = 1.62e2;
+    if (radius >= R) return 0.0;
+    return rho0 * (1.0 - (radius * radius) / (R * R));
+}
+
+double temperature(const mfem::Vector& x) {
+    double radius = x.Norml2();
+    constexpr double R = 5;
+    constexpr double T0 = 1.5e7;
+    if (radius >= R) return 0.0;
+    return T0 * (1.0 - (radius * radius) / (R * R));
+}
+double H1MassFraction(const mfem::Vector &x) {
+    double radius = x.Norml2();
+    // H1 mass fraction, quadratic decrease from 0.7 at center to 0 at surface also with a spiky profile at
+    // mid-radius to simulate burning
+    constexpr double R = 5;
+    if (radius >= R) return 0.0;
+    double radial_part = 0.7 * (1.0 - (radius * radius) / (R * R));
+    double spike = std::exp(-std::pow((radius - R / 2) / 0.5, 2));
+    return std::max(0.0, radial_part - spike);
+
+}
+
+std::unique_ptr<mfem::Mesh> load_mesh(const std::string& filename) {
+    return std::make_unique<mfem::Mesh>(filename.c_str());
+}
+
+void project_radial_function(mfem::GridFunction& gf, std::function<double(const mfem::Vector& x)> g) {
+    mfem::FunctionCoefficient density_coeff([g](const mfem::Vector& x) {
+        return g(x);
+    });
+
+    gf.ProjectCoefficient(density_coeff);
+}
+
+class DiffusionOperator : public mfem::TimeDependentOperator {
+private:
+    std::unique_ptr<mfem::BilinearForm> M, K;
+    mfem::SparseMatrix M_mat, K_mat;
+    mutable std::unique_ptr<mfem::SparseMatrix> T_mat;
+    mutable double current_dt = 0.0;
+    mfem::CGSolver solver;
+    mfem::GSSmoother prec;
+    mutable mfem::Vector rhs; // Removed 'z' as it wasn't used
+public:
+    DiffusionOperator(mfem::FiniteElementSpace& fes, double D_val)
+        : mfem::TimeDependentOperator(fes.GetTrueVSize()) {
+
+        M = std::make_unique<mfem::BilinearForm>(&fes);
+        M->AddDomainIntegrator(new mfem::MassIntegrator());
+        M->Assemble();
+        M->Finalize();
+        M_mat = M->SpMat();
+
+        mfem::ConstantCoefficient D_coeff(D_val);
+        K = std::make_unique<mfem::BilinearForm>(&fes);
+        K->AddDomainIntegrator(new mfem::DiffusionIntegrator(D_coeff));
+        K->Assemble();
+        K->Finalize();
+        K_mat = K->SpMat();
+
+        prec.SetOperator(M_mat);
+        solver.SetOperator(M_mat);
+        solver.SetPreconditioner(prec);
+        solver.SetPrintLevel(1);
+        solver.SetRelTol(1e-8);
+        solver.SetAbsTol(1e-9);
+
+        rhs.SetSize(fes.GetTrueVSize());
+    }
+
+    void Mult(const mfem::Vector& u, mfem::Vector& du_dt) const override {
+        K_mat.Mult(u, rhs);
+        rhs *= -1.0;
+        solver.Mult(rhs, du_dt);
+    }
+
+
+
+    void ImplicitSolve(const double dt, const mfem::Vector &u, mfem::Vector &du_dt) override {
+        if (std::abs(dt - current_dt) > 1e-12 * dt || !T_mat) {
+            T_mat.reset(mfem::Add(1.0, M->SpMat(), dt, K->SpMat()));
+            current_dt = dt;
+            solver.SetOperator(*T_mat);
+            prec.SetOperator(*T_mat);
+        }
+
+        K->SpMat().Mult(u, rhs);
+        rhs *= -1.0;
+
+        solver.Mult(rhs, du_dt);
+    }
+};
+
+
+int main(const int argc, char** argv) {
+    CLI::App app("Diffusion Solver");
+
+    fourdst::config::Config<Options> config;
+    fourdst::config::register_as_cli(config, app);
+    app.parse(argc, argv);
+
+    auto mesh = load_mesh(config->input_mesh);
+
+    const mfem::H1_FECollection fec(1, mesh->Dimension());
+    mfem::FiniteElementSpace fes(mesh.get(), &fec);
+    mfem::GridFunction rho_gf(&fes);
+    mfem::GridFunction temp_gf(&fes);
+    mfem::GridFunction h1_gf(&fes);
+
+
+    project_radial_function(rho_gf, density);
+    project_radial_function(temp_gf, temperature);
+    project_radial_function(h1_gf, H1MassFraction);
+
+
+    double t = 0.0;
+    double t_final = 2e5;
+    double dt = 10;
+    DiffusionOperator op(fes, 1e-3);
+
+    mfem::SDIRK34Solver solver;
+    solver.Init(op);
+
+
+    // ViewMesh(config->vishost, config->visport, *mesh, h1_gf, "Initial H1 Mass Fraction", true);
+
+    mfem::socketstream sol_sock(config->vishost.c_str(), config->visport);
+    size_t step = 0;
+    while (t < t_final) {
+        solver.Step(h1_gf, t, dt);
+        if (step % 10 == 0) {
+            sol_sock << "solution\n" << *mesh << h1_gf;
+            sol_sock << std::flush;
+        }
+        if (step == 0) {
+            sol_sock << "window_title 'Diffusion of H1 Mass Fraction'\n";
+            sol_sock << "keys 'iIzzMaagpmtppc'\n";
+            sol_sock << std::flush;
+            std::cout << "press enter to start the diffusion...";
+            std::cin.get();
+        }
+        step++;
+    }
+    // ViewMesh(config->vishost, config->visport, *mesh, h1_gf, "Final H1 Mass Fraction", true);
+}

meson.build

@@ -0,0 +1,8 @@
+project('libconfig', ['cpp', 'c'], version: 'v2.1.0', default_options: ['cpp_std=c++23'], meson_version: '>=1.5.0')
+
+add_project_arguments('-fvisibility=default', language: 'cpp')
+
+subdir('build-config')
+
+executable('test_diffusion', 'main.cpp', dependencies: [mfem_dep, cli11_dep, config_dep, composition_dep])
+

readme.md

@@ -0,0 +1,22 @@
+# Simple Diffusion Example
+
+Simple example of diffusion using MFEM and stroid. This is built as a test case for me to verify 
+that stroid generates well conditioned meshes and to test implicit time stepping schemes in MFEM.
+
+This model makes many simplifications for the sake of clarity:
+- Constant diffusion coefficient
+- No advection
+- No sources or sinks
+- No cross-diffusion
+- A simple quadratic + ring initial chemical gradient (completely unphysical, just for testing)
+
+## Building
+```bash
+git clone https://github.com/tboudreaux/mfem-diffusion-example.git
+cd mfem-diffusion-example
+meson setup build
+meson compile -C build
+./build/test_diffusion
+```
+
+You will also need glvis running in the background to visualize the output

subprojects/cli11.wrap

@@ -0,0 +1,10 @@
+[wrap-file]
+directory = CLI11-2.6.1
+source_url = https://github.com/CLIUtils/CLI11/archive/refs/tags/v2.6.1.tar.gz
+source_filename = CLI11-2.6.1.tar.gz
+source_hash = 377691f3fac2b340f12a2f79f523c780564578ba3d6eaf5238e9f35895d5ba95
+source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/cli11_2.6.1-1/CLI11-2.6.1.tar.gz
+wrapdb_version = 2.6.1-1
+
+[provide]
+dependency_names = CLI11

subprojects/libcomposition.wrap

@@ -0,0 +1,4 @@
+[wrap-git]
+url = https://github.com/4D-STAR/libcomposition.git
+revision = v2.3.1
+depth = 1

subprojects/mfem.wrap

@@ -0,0 +1,7 @@
+[wrap-git]
+url = https://github.com/mfem/mfem.git
+revision = v4.8-rc0
+diff_files = mfem/disable_mfem_selfcheck.patch
+depth = 1
+
+[cmake]