From 72f5f894d1e65a273e5b1efdc18a08043ecba67b Mon Sep 17 00:00:00 2001 From: cmdoug <68232338+cmdoug@users.noreply.github.com> Date: Mon, 11 May 2026 23:36:22 -0400 Subject: [PATCH 1/7] Add examples for ff-bifbox paper --- examples/README.md | 33 +-- .../README_brusselator.md | 241 ++++++++++++++++++ .../douglas_jolivet_2026/README_cylinder.md | 90 +++++++ .../douglas_jolivet_2026/README_structure.md | 181 +++++++++++++ .../douglas_jolivet_2026/brusselator/.gitkeep | 0 .../douglas_jolivet_2026/cylinder/.gitkeep | 0 .../eqns_douglas_jolivet_2026_brusselator.idp | 84 ++++++ .../eqns_douglas_jolivet_2026_cylinder.idp | 121 +++++++++ .../eqns_douglas_jolivet_2026_structure.idp | 132 ++++++++++ .../extend_brusselator.md | 136 ++++++++++ .../douglas_jolivet_2026/extend_cylinder.md | 59 +++++ .../douglas_jolivet_2026/extend_structure.md | 76 ++++++ .../douglas_jolivet_2026/mesh_brusselator.md | 86 +++++++ .../douglas_jolivet_2026/mesh_cylinder.md | 37 +++ .../douglas_jolivet_2026/mesh_structure.md | 67 +++++ ...tings_douglas_jolivet_2026_brusselator.idp | 59 +++++ ...settings_douglas_jolivet_2026_cylinder.idp | 60 +++++ ...ettings_douglas_jolivet_2026_structure.idp | 60 +++++ .../douglas_jolivet_2026/structure/.gitkeep | 0 19 files changed, 1506 insertions(+), 16 deletions(-) create mode 100644 examples/douglas_jolivet_2026/README_brusselator.md create mode 100644 examples/douglas_jolivet_2026/README_cylinder.md create mode 100644 examples/douglas_jolivet_2026/README_structure.md create mode 100644 examples/douglas_jolivet_2026/brusselator/.gitkeep create mode 100644 examples/douglas_jolivet_2026/cylinder/.gitkeep create mode 100644 examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_brusselator.idp create mode 100644 examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_cylinder.idp create mode 100644 examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_structure.idp create mode 100644 examples/douglas_jolivet_2026/extend_brusselator.md create mode 100644 examples/douglas_jolivet_2026/extend_cylinder.md create mode 100644 examples/douglas_jolivet_2026/extend_structure.md create mode 100644 examples/douglas_jolivet_2026/mesh_brusselator.md create mode 100644 examples/douglas_jolivet_2026/mesh_cylinder.md create mode 100644 examples/douglas_jolivet_2026/mesh_structure.md create mode 100644 examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_brusselator.idp create mode 100644 examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_cylinder.idp create mode 100644 examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp create mode 100644 examples/douglas_jolivet_2026/structure/.gitkeep diff --git a/examples/README.md b/examples/README.md index c0c8aa8..5041e8c 100644 --- a/examples/README.md +++ b/examples/README.md @@ -21,19 +21,20 @@ This folder contains example implementations that reproduce selected results fro | Folder | Description | Citation | | :--- | :--- | :--- | -| [brokof_etal_2024](https://github.com/cmdoug/ff-bifbox/tree/main/examples/brokof_etal_2024) | stability and resolvent analysis of a reacting compressible flow in a duct with acoustic characteristic boundary conditions | Brokof, P., Douglas, C. M., & Polifke, W. (2024). The role of hydrodynamic shear in the thermoacoustic response of slit flames. _Proceedings of the Combustion Institute_, 40(1), 105362. doi:[10.1016/j.proci.2024.105362](https://doi.org/10.1016/j.proci.2024.105362). | -| [chevalier_etal_2024](https://github.com/cmdoug/ff-bifbox/tree/main/examples/chevalier_etal_2024)| resolvent analysis of swirling jet mean flow with Spalart--Allmaras turbulence model | Chevalier, Q., Douglas, C., & Lesshafft, L. (2024). Resolvent analysis of swirling turbulent jets. _Theoretical and Computational Fluid Dynamics_. doi:[10.1007/s00162-024-00704-2](https://doi.org/10.1007/s00162-024-00704-2). | -| [douglas_2024](https://github.com/cmdoug/ff-bifbox/tree/main/examples/douglas_2024) | analysis of various outflow boundary conditions for swirling flows | Douglas, C. M. (2024). A balanced outflow boundary condition for swirling flows. _Theoretical and Compuational Fluid Dynamics_. doi:[10.1007/s00162-024-00701-5](https://doi.org/10.1007/s00162-024-00701-5). | -| [douglas_etal_2021](https://github.com/cmdoug/ff-bifbox/tree/main/examples/douglas_etal_2021) | bifurcation analysis of an incompressible axisymmetric swirling jet with axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2021). Nonlinear dynamics of fully developed swirling jets. _Journal of Fluid Mechanics_, 924, A14. doi:[10.1017/jfm.2021.615](https://doi.org/10.1017/jfm.2021.615). | -| [douglas_etal_2022](https://github.com/cmdoug/ff-bifbox/tree/main/examples/douglas_etal_2022) | bifurcation analysis of an incompressible axisymmetric swirling annular jet with coordinate transformation and axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2022). Dynamics and bifurcations of laminar annular swirling and non-swirling jets. _Journal of Fluid Mechanics_, 943, A35. doi:[10.1017/jfm.2022.453](https://doi.org/10.1017/jfm.2022.453). | -| [douglas_etal_2023](https://github.com/cmdoug/ff-bifbox/tree/main/examples/douglas_etal_2023) | bifurcation analysis of a reacting weakly-compressible axisymmetric jet with axisymmetry breaking | Douglas, C. M., Polifke, W., & Lesshafft, L. (2023). Flash-back, blow-off, and symmetry breaking of premixed conical flames. _Combustion and Flame_, 258(2), 113060. doi:[10.1016/j.combustflame.2023.113060](https://doi.org/10.1016/j.combustflame.2023.113060). | -| [fani_etal_2018](https://github.com/cmdoug/ff-bifbox/tree/main/examples/fani_etal_2018) | bifurcation analysis of a compressible 2-D flow over a cylinder with reflective symmetry breaking | Fani, A., Citro, V., Giannetti, F., & Auteri, F. (2018). Computation of the bluff-body sound generation by a self-consistent mean flow formulation. _Physics of Fluids_, 30(3), 036102. doi:[10.1063/1.4997536](https://doi.org/10.1063/1.4997536). | -| [FK_problem](https://github.com/cmdoug/ff-bifbox/tree/main/examples/FK_problem) | bifurcation analysis of Frank–Kamenetskii thermal runaway (explosion) in a cylinderical vessel | I. B. Zeldovich, G. I. Barenblatt, V. B. Librovich, G. M. Makhviladze, _Mathematical theory of combustion and explosions_, 1985. | -| [garnaud_2012](https://github.com/cmdoug/ff-bifbox/tree/main/examples/garnaud_2012) | resolvent analysis of an incompressible axisymmetric laminar jet | Garnaud, X. (2012). _Modes, transient dynamics and forced response of circular jets_. [Doctoral dissertation, Ecole Polytechnique], HAL ID:[tel-00740133](https://theses.hal.science/tel-00740133). | -| [marquet_larsson_2015](https://github.com/cmdoug/ff-bifbox/tree/main/examples/marquet_larsson_2015) | bifurcation analysis of an incompressible 3-D flow over a thin flat plate with coordinate transformation and multiple symmetries | Marquet, O., & Larsson, M. (2015). Global wake instabilities of low aspect-ratio flat-plates. _European Journal of Mechanics - B/Fluids_, 49(B), 400-412. doi:[10.1016/j.euromechflu.2014.05.005](https://doi.org/10.1016/j.euromechflu.2014.05.005). | -| [meliga_etal_2012](https://github.com/cmdoug/ff-bifbox/tree/main/examples/meliga_etal_2012) | bifurcation analysis of an axisymmetric incompressible Grabowski--Berger vortex with axisymmetry breaking and weakly nonlinear analysis | Meliga, P., Gallaire, F., & Chomaz, J.-M. (2012). A weakly nonlinear mechanism for mode selection in swirling jets. _Journal of Fluid Mechanics_, 699, 216–262. doi:[10.1017/jfm.2012.93](https://doi.org/10.1017/jfm.2012.93). | -| [moreno-boza_etal_2018](https://github.com/cmdoug/ff-bifbox/tree/main/examples/moreno-boza_etal_2018) | analysis of puffing instability in an axisymmetric non-premixed pool flame | Moreno-Boza, D., Coenen, W., Carpio, J., Sánchez A.L., & Williams, F.A. (2018). On the critical conditions for pool-fire puffing. _Combustion and Flame_, 192, 426-438. doi:[10.1016/j.combustflame.2018.02.011](https://doi.org/10.1016/j.combustflame.2018.02.011). | -| [moulin_etal_2019](https://github.com/cmdoug/ff-bifbox/tree/main/examples/moulin_etal_2019) | stability analysis of a 3-D incompressible flow over a thin flat plate using the mAL preconditioner | Moulin, J., Jolivet, P., & Marquet, O. (2019). Augmented Lagrangian preconditioner for large-scale hydrodynamic stability analysis. _Computer Methods in Applied Mechanics and Engineering_, 351, 718-743. doi:[10.1016/j.cma.2019.03.052](https://doi.org/10.1016/j.cma.2019.03.052). | -| [pralits_etal_2010](https://github.com/cmdoug/ff-bifbox/tree/main/examples/pralits_etal_2010) | bifurcation analysis of an incompressible 2-D flow over a rotating cylinder | Pralits, J. O., Brandt, L., & Giannetti, F. (2010). Instability and sensitivity of the flow around a rotating circular cylinder. _Journal of Fluid Mechanics_, 650, 513–536. doi:[10.1017/S0022112009993764](https://doi.org/10.1017/S0022112009993764). | -| [sipp_lebedev_2007](https://github.com/cmdoug/ff-bifbox/tree/main/examples/sipp_lebedev_2007) | bifurcation analysis of the incompressible 2-D flows over a cylinder and open cavity with weakly nonlinear analysis | Sipp, D., & Lebedev, A. (2007). Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. _Journal of Fluid Mechanics_, 593, 333–358. doi:[10.1017/S0022112007008907](https://doi.org/10.1017/S0022112007008907). | -| [wang_etal_2024](https://github.com/cmdoug/ff-bifbox/tree/main/examples/wang_etal_2024) | global linear and nonlinear analysis of an axisymmetric lean premixed V-flame in a weakly-compressible annular jet | Wang, C., Douglas, C., Guan, Y., Xu. C, & Lesshafft, L. (2024). Onset of global instability in a premixed annular V-flame. _Journal of Fluid Mechanics_, 998, A23. doi:[10.1017/jfm.2024.869](https://doi.org/10.1017/jfm.2024.869). | \ No newline at end of file +| [brokof_etal_2024](./brokof_etal_2024) | stability and resolvent analysis of a reacting compressible flow in a duct with acoustic characteristic boundary conditions | Brokof, P., Douglas, C. M., & Polifke, W. (2024). The role of hydrodynamic shear in the thermoacoustic response of slit flames. _Proceedings of the Combustion Institute_, 40(1), 105362. doi:[10.1016/j.proci.2024.105362](https://doi.org/10.1016/j.proci.2024.105362). | +| [chevalier_etal_2024](./chevalier_etal_2024)| resolvent analysis of swirling jet mean flow with Spalart--Allmaras turbulence model | Chevalier, Q., Douglas, C., & Lesshafft, L. (2024). Resolvent analysis of swirling turbulent jets. _Theoretical and Computational Fluid Dynamics_. doi:[10.1007/s00162-024-00704-2](https://doi.org/10.1007/s00162-024-00704-2). | +| [douglas_2024](./douglas_2024) | analysis of various outflow boundary conditions for swirling flows | Douglas, C. M. (2024). A balanced outflow boundary condition for swirling flows. _Theoretical and Computational Fluid Dynamics_. doi:[10.1007/s00162-024-00701-5](https://doi.org/10.1007/s00162-024-00701-5). | +| [douglas_etal_2021](./douglas_etal_2021) | bifurcation analysis of an incompressible axisymmetric swirling jet with axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2021). Nonlinear dynamics of fully developed swirling jets. _Journal of Fluid Mechanics_, 924, A14. doi:[10.1017/jfm.2021.615](https://doi.org/10.1017/jfm.2021.615). | +| [douglas_etal_2022](./douglas_etal_2022) | bifurcation analysis of an incompressible axisymmetric swirling annular jet with coordinate transformation and axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2022). Dynamics and bifurcations of laminar annular swirling and non-swirling jets. _Journal of Fluid Mechanics_, 943, A35. doi:[10.1017/jfm.2022.453](https://doi.org/10.1017/jfm.2022.453). | +| [douglas_etal_2023](./douglas_etal_2023) | bifurcation analysis of a reacting weakly-compressible axisymmetric jet with axisymmetry breaking | Douglas, C. M., Polifke, W., & Lesshafft, L. (2023). Flash-back, blow-off, and symmetry breaking of premixed conical flames. _Combustion and Flame_, 258(2), 113060. doi:[10.1016/j.combustflame.2023.113060](https://doi.org/10.1016/j.combustflame.2023.113060). | +| [douglas_jolivet_2026](./douglas_jolivet_2026) | bifurcation analysis of a 3-D Brusselator system, a 3-D plate buckling system, and a 2-D compressible Navier--Stokes system | Douglas, C. M., & Jolivet, P. (2026). ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs. doi:[10.48550/arXiv.2509.18429](https://doi.org/10.48550/arXiv.2509.18429). | +| [fani_etal_2018](./fani_etal_2018) | bifurcation analysis of a compressible 2-D flow over a cylinder with reflective symmetry breaking | Fani, A., Citro, V., Giannetti, F., & Auteri, F. (2018). Computation of the bluff-body sound generation by a self-consistent mean flow formulation. _Physics of Fluids_, 30(3), 036102. doi:[10.1063/1.4997536](https://doi.org/10.1063/1.4997536). | +| [FK_problem](./FK_problem) | bifurcation analysis of Frank–Kamenetskii thermal runaway (explosion) in a cylinderical vessel | \Zeldovich, I. B., Barenblatt, G. I., Librovich, V. B., & Makhviladze, G. M. _Mathematical theory of combustion and explosions_, 1985. | +| [garnaud_2012](./garnaud_2012) | resolvent analysis of an incompressible axisymmetric laminar jet | Garnaud, X. (2012). _Modes, transient dynamics and forced response of circular jets_. [Doctoral dissertation, Ecole Polytechnique], HAL ID:[tel-00740133](https://theses.hal.science/tel-00740133). | +| [marquet_larsson_2015](./marquet_larsson_2015) | bifurcation analysis of an incompressible 3-D flow over a thin flat plate with coordinate transformation and multiple symmetries | Marquet, O., & Larsson, M. (2015). Global wake instabilities of low aspect-ratio flat-plates. _European Journal of Mechanics - B/Fluids_, 49(B), 400-412. doi:[10.1016/j.euromechflu.2014.05.005](https://doi.org/10.1016/j.euromechflu.2014.05.005). | +| [meliga_etal_2012](./meliga_etal_2012) | bifurcation analysis of an axisymmetric incompressible Grabowski--Berger vortex with axisymmetry breaking and weakly nonlinear analysis | Meliga, P., Gallaire, F., & Chomaz, J.-M. (2012). A weakly nonlinear mechanism for mode selection in swirling jets. _Journal of Fluid Mechanics_, 699, 216–262. doi:[10.1017/jfm.2012.93](https://doi.org/10.1017/jfm.2012.93). | +| [moreno-boza_etal_2018](./moreno-boza_etal_2018) | analysis of puffing instability in an axisymmetric non-premixed pool flame | Moreno-Boza, D., Coenen, W., Carpio, J., Sánchez, A. L., & Williams, F.A. (2018). On the critical conditions for pool-fire puffing. _Combustion and Flame_, 192, 426-438. doi:[10.1016/j.combustflame.2018.02.011](https://doi.org/10.1016/j.combustflame.2018.02.011). | +| [moulin_etal_2019](./moulin_etal_2019) | stability analysis of a 3-D incompressible flow over a thin flat plate using the mAL preconditioner | Moulin, J., Jolivet, P., & Marquet, O. (2019). Augmented Lagrangian preconditioner for large-scale hydrodynamic stability analysis. _Computer Methods in Applied Mechanics and Engineering_, 351, 718-743. doi:[10.1016/j.cma.2019.03.052](https://doi.org/10.1016/j.cma.2019.03.052). | +| [pralits_etal_2010](./pralits_etal_2010) | bifurcation analysis of an incompressible 2-D flow over a rotating cylinder | Pralits, J. O., Brandt, L., & Giannetti, F. (2010). Instability and sensitivity of the flow around a rotating circular cylinder. _Journal of Fluid Mechanics_, 650, 513–536. doi:[10.1017/S0022112009993764](https://doi.org/10.1017/S0022112009993764). | +| [sipp_lebedev_2007](./sipp_lebedev_2007) | bifurcation analysis of the incompressible 2-D flows over a cylinder and open cavity with weakly nonlinear analysis | Sipp, D., & Lebedev, A. (2007). Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. _Journal of Fluid Mechanics_, 593, 333–358. doi:[10.1017/S0022112007008907](https://doi.org/10.1017/S0022112007008907). | +| [wang_etal_2024](./wang_etal_2024) | global linear and nonlinear analysis of an axisymmetric lean premixed V-flame in a weakly-compressible annular jet | Wang, C., Douglas, C. M., Guan, Y., Xu C., & Lesshafft, L. (2024). Onset of global instability in a premixed annular V-flame. _Journal of Fluid Mechanics_, 998, A23. doi:[10.1017/jfm.2024.869](https://doi.org/10.1017/jfm.2024.869). | \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/README_brusselator.md b/examples/douglas_jolivet_2026/README_brusselator.md new file mode 100644 index 0000000..6a8089e --- /dev/null +++ b/examples/douglas_jolivet_2026/README_brusselator.md @@ -0,0 +1,241 @@ +# Examples: Douglas & Jolivet (202x) +This file shows an example `ff-bifbox` workflow for reproducing the results of the study: +```tex +@article{douglas_jolivet_2026, + title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, + author={Douglas, Christopher M. and Jolivet, Pierre R.}, + year={2026}, +} +``` +The commands below illustrate how to perform a bifurcation analysis of the Brusselator in 3-D using `ff-bifbox`. + +## Setup environment for `ff-bifbox` +1. Navigate to the main `ff-bifbox` directory. +```sh +cd ~/your/path/to/ff-bifbox/ +``` +2. Export working directory and number of processors for easy reference. +```sh +export workdir=examples/douglas_jolivet_2026/brusselator +export nproc=4 +``` + +# Example 1: 3-D Brusselator + +1. Create symbolic links for governing equations and solver settings. +```sh +ln -sf examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_brusselator.idp eqns.idp +ln -sf examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_brusselator.idp settings.idp +``` + +2. Build initial meshes +```sh +FreeFem++-mpi -v 0 examples/douglas_jolivet_2026/mesh_brusselator.md -mo $workdir/cube +``` + +3. Compute leading modes over the base state at $L=2$. +```sh + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo sss -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 0,0,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo nss -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 1,0,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo sns -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 0,1,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo nns -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 1,1,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo snn -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 0,1,1 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -mi cube_eighth.mesh -fo nnn -A 2 -B 5.45 -Dx 0.008 -Dy 0.004 -1/L^2 0.25 -eps_target 0.3+2.1i -eps_nev 6 -eps_pos_gen_non_hermitian -sym 1,1,1 +``` + +4. Compute Hopf points +- 1-D solutions: +```sh +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sss.mode -fo brusselator100 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nss.mode -fo brusselator200 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sss_3.mode -fo brusselator300 -param 1/L^2 -snes_rtol 0 +``` +- 2-D solutions: +```sh +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns.mode -fo brusselator110 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sss_1.mode -fo brusselator120 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns_3.mode -fo brusselator130 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nns.mode -fo brusselator210 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nss_1.mode -fo brusselator220 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nns_2.mode -fo brusselator230 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns_2.mode -fo brusselator310 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sss_5.mode -fo brusselator320 -param 1/L^2 -snes_rtol 0 +``` +- 3-D solutions: +```sh +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi snn.mode -fo brusselator111 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nnn.mode -fo brusselator211 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns_1.mode -fo brusselator112 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nns_1.mode -fo brusselator212 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sss_4.mode -fo brusselator122 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi snn_3.mode -fo brusselator311 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi snn_2.mode -fo brusselator131 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nss_3.mode -fo brusselator222 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns_4.mode -fo brusselator312 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi sns_5.mode -fo brusselator132 -param 1/L^2 -snes_rtol 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi nnn_2.mode -fo brusselator231 -param 1/L^2 -snes_rtol 0 +``` + +5. Continue Hopf bifurcations along branches of periodic orbits up to $L=2$ +- 1-D solutions: +```sh +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator100.hopf -fo brusselator100 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator100_2.porb -fo brusselator100 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.25 -dmax 1e10 -mono 0 -contorder 0 -count 2 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator200.hopf -fo brusselator200 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator200_2.porb -fo brusselator200 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.1 -dmax 1e10 -mono 0 -contorder 0 -count 2 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator300.hopf -fo brusselator300 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator300_2.porb -fo brusselator300 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 +``` +- 2-D solutions +```sh +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator110.hopf -fo brusselator110 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator110_2.porb -fo brusselator110 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.2 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator120.hopf -fo brusselator120 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator120_2.porb -fo brusselator120 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator130.hopf -fo brusselator130 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator130_2.porb -fo brusselator130 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.01 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator210.hopf -fo brusselator210 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator210_2.porb -fo brusselator210 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.1 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator220.hopf -fo brusselator220 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -amax 90 -kmax 1e10 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator220_2.porb -fo brusselator220 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator230.hopf -fo brusselator230 -param 1/L^2 -maxcount 2 -h0 -3 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator230_2.porb -fo brusselator230 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.01 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator310.hopf -fo brusselator310 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator310_2.porb -fo brusselator310 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator320.hopf -fo brusselator320 -param 1/L^2 -maxcount 2 -h0 -1 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator320_2.porb -fo brusselator320 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.0025 -dmax 1e10 -mono 0 -contorder 0 -count 2 +``` +- 3-D solutions +```sh +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator111.hopf -fo brusselator111 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator111_2.porb -fo brusselator111 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.1 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator211.hopf -fo brusselator211 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator211_2.porb -fo brusselator211 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.1 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator112.hopf -fo brusselator112 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator112_2.porb -fo brusselator112 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator212.hopf -fo brusselator212 -param 1/L^2 -maxcount 2 -h0 -5 -Nh 6 -amax 90 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator212_2.porb -fo brusselator212 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.05 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator122.hopf -fo brusselator122 -param 1/L^2 -maxcount 1 -h0 -5 -Nh 6 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator122_1.porb -fo brusselator122 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.02 -dmax 1e10 -mono 0 -contorder 0 -count 1 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator311.hopf -fo brusselator311 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator311_2.porb -fo brusselator311 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.02 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator131.hopf -fo brusselator131 -param 1/L^2 -maxcount 2 -h0 -2 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator131_2.porb -fo brusselator131 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.01 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator222.hopf -fo brusselator222 -param 1/L^2 -maxcount 2 -h0 -1 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator222_2.porb -fo brusselator222 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.005 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator312.hopf -fo brusselator312 -param 1/L^2 -maxcount 2 -h0 -1 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator312_2.porb -fo brusselator312 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.00125 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator132.hopf -fo brusselator132 -param 1/L^2 -maxcount 2 -h0 -1 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator132_2.porb -fo brusselator132 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.0025 -dmax 1e10 -mono 0 -contorder 0 -count 2 + +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator231.hopf -fo brusselator231 -param 1/L^2 -maxcount 2 -h0 -1 -Nh 6 -kmax 1e10 -amax 90 -mono 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator231_2.porb -fo brusselator231 -param 1/L^2 -paramtarget 0.25 -maxcount -1 -h0 -0.0025 -dmax 1e10 -mono 0 -contorder 0 -count 2 +``` + +6. (OPTIONAL) Floquet analysis calculations along the 1-D solution branches +```sh +cd "$workdir" && declare -a porblist=(brusselator100_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 1,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100 -eps_target 0.3+2.1i -sym 1,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done + +cd "$workdir" && declare -a porblist=(brusselator200_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 1,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200 -eps_target 0.3+2.1i -sym 1,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done + +cd "$workdir" && declare -a porblist=(brusselator300_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 1,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300 -eps_target 0.3+2.1i -sym 1,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done +``` + +6. Compute other 1-D branches using deflation and continue them +```sh +ff-mpirun -np $nproc porbcompute.md -v 0 -dir $workdir -fi brusselator100_10.porb -fo brusselator100_nss -1/L^2 0.65 +ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi brusselator100_nss.porb -fo brusselator100_nss -eps_target 0.01+1.81i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 5 +ff-mpirun examples/douglas_jolivet_2026/extend_brusselator.md -v 0 -dir $workdir -fi brusselator100_nss.floq -fo brusselator100_nss -amp 0.1 +ff-mpirun -np $nproc porbdeflate.md -v 0 -dir $workdir -mi cube_Qx.mesh -fi brusselator100_nssguess.porb -fi2 brusselator100_nssstart.porb -fo brusselator100_nss -ndeflate 1 -snes_rtol 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator100_nss-1.porb -fo brusselator100_nss -param 1/L^2 -h0 1 -paramtarget 0.25 -maxcount -1 + +ff-mpirun -np $nproc porbcompute.md -v 0 -dir $workdir -fi brusselator100_12.porb -fo brusselator100_nss2 -1/L^2 0.28 +ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi brusselator100_nss2.porb -fo brusselator100_nss2 -eps_target 0.01+1.84i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 5 +ff-mpirun examples/douglas_jolivet_2026/extend_brusselator.md -v 0 -dir $workdir -fi brusselator100_nss2.floq -fo brusselator100_nss2 -amp 0.1 +ff-mpirun -np $nproc porbdeflate.md -v 0 -dir $workdir -mi cube_Qx.mesh -fi brusselator100_nss2guess.porb -fi2 brusselator100_nss2start.porb -fo brusselator100_nss2 -ndeflate 1 -snes_rtol 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator100_nss2-1.porb -fo brusselator100_nss2 -param 1/L^2 -h0 -1 -paramtarget 0.25 -maxcount -1 + +ff-mpirun -np $nproc porbcompute.md -v 0 -dir $workdir -fi brusselator200_5.porb -fo brusselator200_sss -1/L^2 0.57 +ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi brusselator200_sss.porb -fo brusselator200_sss -eps_target 0+1.98i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 5 +ff-mpirun examples/douglas_jolivet_2026/extend_brusselator.md -v 0 -dir $workdir -fi brusselator200_sss_2.floq -fo brusselator200_sss -amp 0.1 +ff-mpirun -np $nproc porbdeflate.md -v 0 -dir $workdir -mi cube_Qx.mesh -fi brusselator200_sssguess.porb -fi2 brusselator200_sssstart.porb -fo brusselator200_sss -ndeflate 1 -snes_rtol 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator200_sss-1.porb -fo brusselator200_sss -param 1/L^2 -h0 -0.005 -paramtarget 0.25 -maxcount -1 -contorder 0 -mono 0 -dmax 1e10 + +ff-mpirun -np $nproc porbcompute.md -v 0 -dir $workdir -fi brusselator300_4.porb -fo brusselator300_nss -1/L^2 0.285 +ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi brusselator300_nss.porb -fo brusselator300_nss -eps_target 0+2.01i -sym 1,0,0 -eps_pos_gen_non_hermitian -eps_nev 5 +ff-mpirun examples/douglas_jolivet_2026/extend_brusselator.md -v 0 -dir $workdir -fi brusselator300_nss_3.floq -fo brusselator300_nss -amp 0.1 +ff-mpirun -np $nproc porbdeflate.md -v 0 -dir $workdir -mi cube_Qx.mesh -fi brusselator300_nssguess.porb -fi2 brusselator300_nssstart.porb -fo brusselator300_nss -ndeflate 1 -snes_rtol 0 +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi brusselator300_nss-1.porb -fo brusselator300_nss -param 1/L^2 -h0 -0.005 -paramtarget 0.25 -maxcount -1 -contorder 0 -mono 0 -dmax 1e10 +``` + +6. (OPTIONAL) Floquet analysis calculations along the 1-D solution branches +```sh + +cd "$workdir" && declare -a porblist=(brusselator100_nss_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss -eps_target 0.3+1.9i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss -eps_target 0.3+1.9i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss -eps_target 0.3+1.9i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done + +cd "$workdir" && declare -a porblist=(brusselator100_nss2_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss2 -eps_target 0.3+1.9i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss2 -eps_target 0.3+1.9i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator100_nss2 -eps_target 0.3+1.9i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done + +cd "$workdir" && declare -a porblist=(brusselator200_sss_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200_sss -eps_target 0.3+1.85i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200_sss -eps_target 0.3+1.85i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 10 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator200_sss -eps_target 0.3+1.85i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 10 +done + +cd "$workdir" && declare -a porblist=(brusselator300_nss_*[0-9].porb) && cd - +for porbfile in "${porblist[@]}"; do + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300_nss -eps_target 0.3+2i -sym 0,0,0 -eps_pos_gen_non_hermitian -eps_nev 16 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300_nss -eps_target 0.3+2i -sym 0,1,0 -eps_pos_gen_non_hermitian -eps_nev 16 + ff-mpirun -np $nproc floqcompute.md -v 0 -dir $workdir -fi $porbfile -so brusselator300_nss -eps_target 0.3+2i -sym 0,1,1 -eps_pos_gen_non_hermitian -eps_nev 16 +done +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/README_cylinder.md b/examples/douglas_jolivet_2026/README_cylinder.md new file mode 100644 index 0000000..e154f14 --- /dev/null +++ b/examples/douglas_jolivet_2026/README_cylinder.md @@ -0,0 +1,90 @@ +# Compressible Flow Example: +This file shows an example `ff-bifbox` workflow for reproducing the results of the paper: +```tex +@article{douglas_jolivet_2026, + title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, + author={Douglas, Christopher M. and Jolivet, Pierre R.}, + year={2026}, +} +``` +The commands below illustrate how to analyze a 2-D compressible flow past a cylinder using `ff-bifbox`. + +## Setup environment for `ff-bifbox` +1. Navigate to the main `ff-bifbox` directory. +```sh +cd ~/your/path/to/ff-bifbox/ +``` +2. Export working directory and number of processors for easy reference. +```sh +export workdir=examples/douglas_jolivet_2026/cylinder_alt +export nproc=4 +``` + +# Example 3: compressible flow past a cylinder + +1. Create symbolic links for governing equations and solver settings. +```sh +ln -sf examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_cylinder_alt.idp eqns.idp +ln -sf examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_cylinder_alt.idp settings.idp +``` + +2. Build initial mesh +```sh +FreeFem++-mpi -v 0 examples/douglas_jolivet_2026/mesh_cylinder.md -mo $workdir/cylinder +``` + +3. Compute base states on the created meshes at $Re=40,50,70,90$ and $Ma=0,0.4,0.6,0.8$ from default guess +```sh +ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -mi cylinder.msh -fo cylinder_Ma0p0_Re40 -1/Re 0.025 -1/Pr 1.38888888889 -Ma^2 0 -gamma 1.4 +ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re40.base -fo cylinder_Ma0p0_Re50 -1/Re 0.02 +ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re50.base -fo cylinder_Ma0p0_Re70 -1/Re 0.014285714285714285 +ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re70.base -fo cylinder_Ma0p0_Re90 -1/Re 0.0111111111111111111 +for Re in 40 50 70 90; do + ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re"$Re".base -fo cylinder_Ma0p4_Re"$Re" -Ma^2 0.16 + ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p4_Re"$Re".base -fo cylinder_Ma0p6_Re"$Re" -Ma^2 0.36 + ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p6_Re"$Re".base -fo cylinder_Ma0p8_Re"$Re" -Ma^2 0.64 + for Ma in 0 4 6 8; do + ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma"_Re"$Re".base -fo cylinder_Ma0p"$Ma"_Re"$Re" -mo cylinder_Ma0p"$Ma"_Re"$Re" -thetamax 0.1 -hmin 1.e-5 -hmax 2 + done +done +``` + +4. Continue base states along $Re$ or $Ma$ with adaptive remeshing +```sh +for Ma in 0 4 6 8; do + ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma"_Re40.base -fo cylinder_Ma0p"$Ma"_Re20 -1/Re 0.05 + ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma"_Re20.base -fo cylinder_Ma0p"$Ma" -param 1/Re -h0 -1 -scount 2 -maxcount -1 -paramtarget 0.01 -mo cylinder_Ma0p"$Ma" -thetamax 0.1 -hmin 1.e-5 -hmax 2 +done +for Re in 40 50 70 90; do + ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re"$Re".base -fo cylinder_Re"$Re" -param Ma^2 -h0 0.25 -scount 2 -maxcount -1 -paramtarget 0.81 -mo cylinder_Re"$Re" -thetamax 0.1 -hmin 1.e-5 -hmax 2 +done +``` + +5. Compute Hopf bifurcation at $Ma=0$ and trace it along the neutral curve +```sh +ff-mpirun -np $nproc basecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re50.base -fo cylinder_Ma0p0_Re47p6 -1/Re 0.021 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re47p6.base -fo cylinder_Ma0p0_Re47p6 -eps_target 0.1+0.7i -sym 1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0_Re47p6.mode -fo cylinder_Ma0p0 -param 1/Re -nf 0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_Ma0p0.hopf -fo cylinder_Ma0p0 -mo cylinder_Ma0p0 -param 1/Re -thetamax 0.1 -hmin 1e-5 -hmax 2 +ff-mpirun -np $nproc hopfcontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p0.hopf -fo cylinder -mo cylinder -thetamax 0.1 -hmin 1e-5 -hmax 2 -param 1/Re -param2 Ma^2 -h0 0.1 -scount 3 -paramtarget 0.01 -maxcount -1 +``` + +6. Compute Hopf bifurcations at $Ma=0.4,0.6,0.8$ and $Re=50,70,90$ +```sh +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_12.hopf -fo cylinder_Ma0p4 -param 1/Re -Ma^2 0.16 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_15specialpt.hopf -fo cylinder_Ma0p6 -param 1/Re -Ma^2 0.36 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_21.hopf -fo cylinder_Ma0p8 -param 1/Re -Ma^2 0.64 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_15specialpt.hopf -fo cylinder_Re50 -param Ma^2 -1/Re 0.02 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_18.hopf -fo cylinder_Re70 -param Ma^2 -1/Re 0.014285714285714285 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_21.hopf -fo cylinder_Re90 -param Ma^2 -1/Re 0.0111111111111111111 +``` + +7. Continue the branches of periodic solutions emanating from the Hopf points along $Re$ and $Ma$. +```sh +for Ma in 0 4 6 8; do +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma".hopf -fo cylinder_Ma0p"$Ma" -mo cyl_Ma0p"$Ma" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param 1/Re -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 +done +for Re in 50 70 90; do +ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Re"$Re".hopf -fo cylinder_Re"$Re" -mo cyl_Re"$Re" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param Ma^2 -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 +done +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/README_structure.md b/examples/douglas_jolivet_2026/README_structure.md new file mode 100644 index 0000000..8fb3ad9 --- /dev/null +++ b/examples/douglas_jolivet_2026/README_structure.md @@ -0,0 +1,181 @@ +# Examples: Douglas & Jolivet (202x) +This file shows an example `ff-bifbox` workflow for reproducing the results of the study: +```tex +@article{douglas_jolivet_2026, + title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, + author={Douglas, Christopher M. and Jolivet, Pierre R.}, + year={2026}, +} +``` +The commands below illustrate how to perform a bifurcation analysis of buckling 3-D structure using `ff-bifbox`. + +## Setup environment for `ff-bifbox` +1. Navigate to the main `ff-bifbox` directory. +```sh +cd ~/your/path/to/ff-bifbox/ +``` +2. Export working directory and number of processors for easy reference. +```sh +export workdir=examples/douglas_jolivet_2026/structure_test +export nproc=4 +``` + +# Example 1: 3-D buckling plate + +1. Create symbolic links for governing equations and solver settings. +```sh +ln -sf examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_structure.idp eqns.idp +ln -sf examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp settings.idp +``` + +2. Build meshes +```sh +FreeFem++-mpi -v 0 examples/douglas_jolivet_2026/mesh_structure.md -mo $workdir/structure +``` + +3. Compute $S$ branch at various $P$ values using predictor--corrector continuation method +```sh +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_S.mesh -fo structure_S -param P -nu 0.3 -h0 0.1 -tgv -2 -paramtarget 1.5e-7 -maxcount -1 +``` + +4. Compute fold points on the $S$ curve +```sh +cd "$workdir" && declare -a foldguesslist=(structure_S_*specialpt.base) && cd - +for guess in "${foldguesslist[@]}"; do + out=$(echo "$guess" | awk -F'specialpt' '{print $1}') + ff-mpirun -np $nproc foldcompute.md -v 0 -dir $workdir -fi $guess -fo $out -param P -tgv -2 +done +``` + +5. (OPTIONAL) Compute the relevant portion of the eigenvalue spectrum of the $S$ branch at each point for each possible symmetry (i.e. $S$, $H_x$, $H_y$, $R_z$). +```sh +cd "$workdir" && declare -a baselist=(structure_S_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_SS -eps_target 0.1+0i -eps_nev 3 -eps_gen_hermitian -sym 0,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_SHx -eps_target 0.1+0i -eps_nev 3 -eps_gen_hermitian -sym 0,1 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_SHy -eps_target 0.1+0i -eps_nev 3 -eps_gen_hermitian -sym 1,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_SRz -eps_target 0.1+0i -eps_nev 3 -eps_gen_hermitian -sym 1,1 +done +``` + +6. Compute eigenmodes near exchanges of stability and iterate to find pitchfork bifurcations associated with symmetry breaking of the $S$ state +```sh +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_4.base -fo structure_SHx -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 0,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SHx.mode -fo structure_SHx -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_30.base -fo structure_SHx2 -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 0,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SHx2.mode -fo structure_SHx2 -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_9.base -fo structure_SHy -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SHy.mode -fo structure_SHy -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_20.base -fo structure_SHy2 -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,0 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SHy2.mode -fo structure_SHy2 -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_14.base -fo structure_SRz -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SRz.mode -fo structure_SRz -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_S_22.base -fo structure_SRz2 -eps_target 0.1+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_SRz2.mode -fo structure_SRz2 -param P -tgv -2 -zero 1 +``` + +7. Extend the computed bifurcations across the domain according to their symmetry and trace them along $P$ +```sh +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_SHx.hopf -v 0 -fo structure_SHx +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_Hx.mesh -fi structure_SHxstart.base -fi2 structure_SHxbranch.base -fo structure_Hx -param P -h0 -0.0071825 -tgv -2 -snes_rtol 0 -mono 0 -maxcount 35 +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_SHy.hopf -v 0 -fo structure_SHy +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_Hy.mesh -fi structure_SHystart.base -fi2 structure_SHybranch.base -fo structure_Hy -param P -h0 -0.03125 -tgv -2 -snes_rtol 0 -mono 0 -maxcount 40 +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_SRz.hopf -v 0 -fo structure_SRz +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_Rz.mesh -fi structure_SRzstart.base -fi2 structure_SRzbranch.base -fo structure_Rz -param P -h0 0.015625 -tgv -2 -snes_rtol 0 -maxcount 41 +``` + +8. Compute fold points on the bifurcation curves +```sh +cd "$workdir" && declare -a foldguesslist=(structure_[H,R][x,y,z]_*specialpt.base) && cd - +for guess in "${foldguesslist[@]}"; do + out=$(echo "$guess" | awk -F'specialpt' '{print $1}') + ff-mpirun -np $nproc foldcompute.md -v 0 -dir $workdir -fi $guess -fo $out -param P -tgv -2 +done +ff-mpirun -np $nproc foldcompute.md -v 0 -dir $workdir -fi structure_Hy_34.base -fo structure_Hy_35 -param P -tgv -2 +``` + +9. (OPTIONAL) Compute the relevant portion of the eigenvalue spectrum of the $R_z$ branch at each point for each broken symmetry. +```sh +cd "$workdir" && declare -a baselist=(structure_Hx_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_HxHx -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 0,1 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_HxA -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +cd "$workdir" && declare -a baselist=(structure_Hy_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_HyHy -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_HyA -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +cd "$workdir" && declare -a baselist=(structure_Rz_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_RzRz -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 0,0 + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_RzA -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +``` + +10. Compute eigenmodes near exchanges of stability and iterate to find pitchfork bifurcations associated with symmetry breaking of the $H_x$, $H_y$, and $R_z$ states +- Bifurcations from $H_x$ to $A$ states: +```sh +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hx_14.base -fo structure_HxA -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HxA.mode -fo structure_HxA -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hx_16.base -fo structure_HxA2 -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HxA2.mode -fo structure_HxA2 -param P -tgv -2 -zero 1 +``` +- Bifurcations from $H_y$ to $A$ states: +```sh +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hy_17.base -fo structure_HyA -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HyA.mode -fo structure_HyA -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hy_26.base -fo structure_HyA2 -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HyA2.mode -fo structure_HyA2 -param P -tgv -2 -zero 1 +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hy_34.base -fo structure_HyA3 -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HyA3.mode -fo structure_HyA3 -param P -tgv -2 -zero 1 +``` +- Bifurcations from $R_z$ to $A$ state: +```sh +ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Rz_23.base -fo structure_RzA -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 +ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_RzA.mode -fo structure_RzA -param P -tgv -2 -zero 1 +``` + +11. Trace bifurcated solutions along $P$ +- From $H_x$ to $A$ +```sh +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_HxA.hopf -v 0 -fo structure_HxA -pv 1 +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_A.mesh -fi structure_HxAstart.base -fi2 structure_HxAbranch.base -fo structure_A -param P -h0 -0.015625 -tgv -2 -snes_rtol 0 -mono 0 -maxcount 29 +``` +- From $H_y$ to $A$ +```sh +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_HyA.hopf -v 0 -fo structure_HyA +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_A.mesh -fi structure_HyAstart.base -fi2 structure_HyAbranch.base -fo structure_A2 -param P -h0 0.015625 -tgv -2 -snes_rtol 0 -mono 0 -maxcount 36 +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_HyA3.hopf -v 0 -fo structure_HyA3 -pv 1 +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_A.mesh -fi structure_HyA3start.base -fi2 structure_HyA3branch.base -fo structure_A3 -param P -h0 0.0078125 -tgv -2 -snes_rtol 0 -mono 0 -maxcount 20 +``` +- From $R_z$ to $A$ +```sh +ff-mpirun -np 1 examples/douglas_jolivet_2026/extend_structure.md -dir $workdir -fi structure_RzA.hopf -v 0 -fo structure_RzA +ff-mpirun -np $nproc basecontinue.md -v 0 -dir $workdir -mi structure_A.mesh -fi structure_RzAstart.base -fi2 structure_RzAbranch.base -fo structure_A4 -param P -h0 -0.03125 -tgv -2 -snes_rtol 0 -maxcount 23 +``` + +12. Compute fold points on the bifurcation curves +```sh +cd "$workdir" && declare -a foldguesslist=(structure_A*specialpt.base) && cd - +for guess in "${foldguesslist[@]}"; do + out=$(echo "$guess" | awk -F'specialpt' '{print $1}') + ff-mpirun -np $nproc foldcompute.md -v 0 -dir $workdir -fi $guess -fo $out -param P -tgv -2 +done +``` + +13. (OPTIONAL) Compute the relevant portion of the eigenvalue spectrum +```sh +cd "$workdir" && declare -a baselist=(structure_A_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_A -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +cd "$workdir" && declare -a baselist=(structure_A2_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_A2 -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +cd "$workdir" && declare -a baselist=(structure_A3_*[0-9].base) && cd - +for baseflow in "${baselist[@]}"; do + ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi $baseflow -so structure_A3 -eps_target 0.1+0i -eps_nev 5 -eps_gen_hermitian -sym 1,1 +done +``` diff --git a/examples/douglas_jolivet_2026/brusselator/.gitkeep b/examples/douglas_jolivet_2026/brusselator/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/examples/douglas_jolivet_2026/cylinder/.gitkeep b/examples/douglas_jolivet_2026/cylinder/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_brusselator.idp b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_brusselator.idp new file mode 100644 index 0000000..72b236d --- /dev/null +++ b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_brusselator.idp @@ -0,0 +1,84 @@ +// +// eqns_douglas_jolivet_2026_brusselator.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// macros + macro vdotu(v, u) (v*u + v#Y*u#Y) //EOM + + macro diff(v, u) (dx(v)*dx(u) + dy(v)*dy(u) + dz(v)*dz(u)) //EOM + +// Define KSP parameters + string KSPparams = "-ksp_type preonly -pc_type lu"; + +// Boundary conditions + macro BoundaryConditions(u, U) + on(BCdX, BCnX, (abs(int(sym(0))) % 2)*BCsX, (abs(int(sym(1))) % 2)*BCsY, (abs(int(sym(2))) % 2)*BCsZ, u = U - params["A"], u#Y = U#Y - params["B"]/params["A"]) + // EOM + macro HomBoundaryConditions(u) + on(BCdX, BCnX, (abs(int(sym(0))) % 2)*BCsX, (abs(int(sym(1))) % 2)*BCsY, (abs(int(sym(2))) % 2)*BCsZ, u = 0, u#Y = 0) + // EOM + +// RESIDUAL OPERATOR + varf vR(defu(dum), defu(v)) + = int3d(Th, qforder=3)( + (ub - params["A"])*v + (vY - v)*(ub*ubY - params["B"])*ub + params["1/L^2"]*(params["Dx"]*diff(v, ub) + params["Dy"]*diff(vY, ubY)) + ) + + BoundaryConditions(dum, ub); + +// JACOBIAN OPERATOR + varf vJ(defu(dum), defu(v)) + = int3d(Th, qforder=3)( + iomega*vdotu(v, dum) + dum*v + (vY - v)*((ub*dumY + 2.0*dum*ubY)*ub - params["B"]*dum) + params["1/L^2"]*(params["Dx"]*diff(v, dum) + params["Dy"]*diff(vY, dumY)) + ) + + int3d(Th, qforder=3)( + iomega*vdotu(v, um) + um*v + (vY - v)*((ub*umY + 2.0*um*ubY)*ub - params["B"]*um) + params["1/L^2"]*(params["Dx"]*diff(v, um) + params["Dy"]*diff(vY, umY)) + ) + + HomBoundaryConditions(dum); + +// MASS OPERATOR + varf vM(defu(dum), defu(v)) + = int3d(Th, qforder=3)( vdotu(v, dum) ) + + int3d(Th, qforder=3)( vdotu(v, um) ) + + HomBoundaryConditions(dum); + +// MASS DERIVATIVE OPERATOR + varf vdM(defu(dum), defu(v)) + = HomBoundaryConditions(dum); + +// MASS 2ND DERIVATIVE OPERATOR + varf vddM(defu(dum), defu(v)) + = HomBoundaryConditions(dum); + +// FORCING/RESPONSE OPERATORS + varf vMq(defu(dum), defu(v)) + = int3d(Th, qforder=3)( vdotu(v, dum) ) + + int3d(Th, qforder=3)( vdotu(v, um) ); + + varf vMf(deff(fm), deff(v)) + = int3d(Th, qforder=3)( vdotu(v, fm) ); + + varf vP(deff(fm), defu(v)) + = int3d(Th, qforder=3)( vdotu(v, fm) ); + +// HESSIAN OPERATOR + varf vH(defu(dum), defu(v)) + = int3d(Th, qforder=3)( + 2.0*(vY - v)*(dum*(ub*umY + um*ubY) + ub*um*dumY) + ) + + int3d(Th, qforder=3)( + 2.0*(vY - v)*(um2*(ub*umY + um*ubY) + ub*um*um2Y) + ) + + HomBoundaryConditions(dum); + +// TRESSIAN OPERATOR + varf vT(defu(dum), defu(v)) + = int3d(Th, qforder=3)( + 2.0*(vY - v)*(um2*(dum*umY + um*dumY) + dum*um*um2Y) + ) + + int3d(Th, qforder=3)( + 2.0*(vY - v)*(um2*(um3*umY + um*um3Y) + um3*um*um2Y) + ) + + HomBoundaryConditions(dum); + + varf vQ(defu(dum), defu(v)) = HomBoundaryConditions(dum); \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_cylinder.idp b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_cylinder.idp new file mode 100644 index 0000000..8886840 --- /dev/null +++ b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_cylinder.idp @@ -0,0 +1,121 @@ +// +// eqns_douglas_jolivet_2026_cylinder.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// macros + func bsds = ( 1.0e-4*(max(x - 50., -20. - x, 0.)^2. + max(y - 20., 0.)^2.) ); //distance to physical domain + macro vdotu(v, u) ( v*u + v#y*u#y ) // velocity inner product + macro div(u) ( dx(u) + dy(u#y) ) // velocity divergence + macro ugradu(v, U, u) ( v*(U*dx(u) + U#y*dy(u)) + v#y*(U*dx(u#y) + U#y*dy(u#y)) ) // velocity advection term + macro visc(v, u) ( (dx(v) - dy(v#y))*(dx(u) - dy(u#y)) + (dy(v) + dx(v#y))*(dx(u#y) + dy(u)) ) // viscous term + macro visc2(v, T, u) ( (dx(T)*v - dy(T)*v#y)*(dx(u) - dy(u#y)) + (dy(T)*v + dx(T)*v#y)*(dx(u#y) + dy(u)) ) // other viscous term + macro diff(g, f) ( dx(g)*dx(f) + dy(g)*dy(f) ) // scalar diffusion term + macro ugradf(U, f) ( U*dx(f) + U#y*dy(f) ) // scalar advection term + +// Define KSP parameters + string KSPparams = "-ksp_type preonly -pc_type lu"; + +// Boundary conditions + macro BoundaryConditions(u, U) + on(BCwall, u = U, u#y = U#y) + + on(BCin, u = U - 1.0, u#y = U#y, u#T = U#T - 1.0) + + on(BClat, BCsym, u#y = U#y) + // EOM + macro HomBoundaryConditions(u) + on(BCwall, u = 0, u#y = 0) + + on(BCin, u = 0, u#y = 0, u#T = 0) + + on(BClat, (abs(int(sym(0))) % 2 == 0)*BCsym, u#y = 0) + + on((abs(int(sym(0))) % 2 != 0)*BCsym, u = 0, u#T = 0) + // EOM + +// RESIDUAL OPERATOR + varf vR(defu(dum), defu(v)) + = int2d(Th)( + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*ugradu(v, ub, ub) - (ubT*div(v) + ugradf(v, ubT))*ubp + params["1/Re"]*(ubT*visc(v, ub) + visc2(v, ubT, ub)) + + vT*((params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(ub, ubT) + (params["gamma"] - 1.0)*ubT*div(ub)) - params["gamma"]*params["Ma^2"]*(params["gamma"] - 1.0)*params["1/Re"]*ubT*visc(ub, ub)) + params["gamma"]*params["1/Re"]*params["1/Pr"]*(ubT*diff(vT, ubT) + vT*diff(ubT, ubT)) + + vp*(params["gamma"]*params["Ma^2"]*ubT*ugradf(ub, ubp) - (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(ub, ubT) - ubT*div(ub))) + + bsds*(vdotu(v, ub) - v + params["gamma"]*((ubT - 1.0)*vT + params["Ma^2"]*ubp*vp)) + ) + + BoundaryConditions(dum, ub); + +// JACOBIAN OPERATOR + varf vJ(defu(dum), defu(v)) + = int2d(Th)( + iomega*((params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, um) + (vT - vp)*umT) + params["gamma"]*params["Ma^2"]*ubT*vp*ump) + + params["gamma"]*params["Ma^2"]*ump*ugradu(v, ub, ub) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradu(v, um, ub) + ugradu(v, ub, um)) - (ubT*div(v) + ugradf(v, ubT))*ump - (umT*div(v) + ugradf(v, umT))*ubp + params["1/Re"]*(umT*visc(v, ub) + ubT*visc(v, um) + visc2(v, umT, ub) + visc2(v, ubT, um)) + + vT*(params["gamma"]*params["Ma^2"]*(ump*(ugradf(ub, ubT) + (params["gamma"] - 1.0)*ubT*div(ub)) - (params["gamma"] - 1.0)*params["1/Re"]*(umT*visc(ub, ub) + ubT*(visc(um, ub) + visc(ub, um)))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(um, ubT) + ugradf(ub, umT) + (params["gamma"] - 1.0)*(umT*div(ub) + ubT*div(um)))) + params["gamma"]*params["1/Re"]*params["1/Pr"]*(umT*diff(vT, ubT) + ubT*diff(vT, umT) + 2.0*vT*diff(ubT, umT)) + + vp*(params["gamma"]*params["Ma^2"]*(umT*ugradf(ub, ubp) + ubT*(ugradf(um, ubp) + ugradf(ub, ump)) - ump*(ugradf(ub, ubT) - ubT*div(ub))) - (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(um, ubT) + ugradf(ub, umT) - umT*div(ub) - ubT*div(um))) + + bsds*(vdotu(v, um) + params["gamma"]*(umT*vT + params["Ma^2"]*ump*vp)) + ) + + int2d(Th)( + iomega*((params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, dum) + (vT - vp)*dumT) + params["gamma"]*params["Ma^2"]*ubT*vp*dump) + + params["gamma"]*params["Ma^2"]*dump*ugradu(v, ub, ub) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradu(v, dum, ub) + ugradu(v, ub, dum)) - (ubT*div(v) + ugradf(v, ubT))*dump - (dumT*div(v) + ugradf(v, dumT))*ubp + params["1/Re"]*(dumT*visc(v, ub) + ubT*visc(v, dum) + visc2(v, dumT, ub) + visc2(v, ubT, dum)) + + vT*(params["gamma"]*params["Ma^2"]*(dump*(ugradf(ub, ubT) + (params["gamma"] - 1.0)*ubT*div(ub)) - (params["gamma"] - 1.0)*params["1/Re"]*(dumT*visc(ub, ub) + ubT*(visc(dum, ub) + visc(ub, dum)))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(dum, ubT) + ugradf(ub, dumT) + (params["gamma"] - 1.0)*(dumT*div(ub) + ubT*div(dum)))) + params["gamma"]*params["1/Re"]*params["1/Pr"]*(dumT*diff(vT, ubT) + ubT*diff(vT, dumT) + 2.0*vT*diff(ubT, dumT)) + + vp*(params["gamma"]*params["Ma^2"]*(dumT*ugradf(ub, ubp) + ubT*(ugradf(dum, ubp) + ugradf(ub, dump)) - dump*(ugradf(ub, ubT) - ubT*div(ub))) - (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(dum, ubT) + ugradf(ub, dumT) - dumT*div(ub) - ubT*div(dum))) + + bsds*(vdotu(v, dum) + params["gamma"]*(dumT*vT + params["Ma^2"]*dump*vp)) + ) + + HomBoundaryConditions(dum); + +// MASS OPERATOR + varf vM(defu(dum), defu(v)) + = int2d(Th)( (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, dum) + (vT - vp)*dumT) + params["gamma"]*params["Ma^2"]*ubT*vp*dump ) + + int2d(Th)( (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, um) + (vT - vp)*umT) + params["gamma"]*params["Ma^2"]*ubT*vp*ump ) + + HomBoundaryConditions(dum); + +// MASS DERIVATIVE OPERATOR + varf vdM(defu(dum), defu(v)) + = int2d(Th)( params["gamma"]*params["Ma^2"]*(dump*(vdotu(v, um) + (vT - vp)*umT) + dumT*vp*ump ) ) + + int2d(Th)( params["gamma"]*params["Ma^2"]*(um2p*(vdotu(v, um) + (vT - vp)*umT) + um2T*vp*ump ) ) + + HomBoundaryConditions(dum); + +// MASS 2ND DERIVATIVE OPERATOR + varf vddM(defu(dum), defu(v)) + = HomBoundaryConditions(dum); + +// FORCING/RESPONSE OPERATORS + varf vMq(defu(dum), defu(v)) + = int2d(Th)( + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, dum) + vp*dumT) + params["Ma^2"]*(vT - params["gamma"]*vp)*dump + ) + + int2d(Th)( + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(vdotu(v, um) + vp*umT) + params["Ma^2"]*(vT - params["gamma"]*vp)*ump + ) + + HomBoundaryConditions(dum); + + varf vMf(deff(fm), deff(v)) + = int2d(Th)( vdotu(v, fm) ); + + varf vP(deff(fm), defu(v)) + = int2d(Th)( vdotu(v, fm) ); + +// HESSIAN OPERATOR + varf vH(defu(dum), defu(v)) + = int2d(Th)( + params["gamma"]*params["Ma^2"]*(iomega*(um2p*(vdotu(v, um) + (vT - vp)*umT) + um2T*vp*ump ) + iomega2*(ump*(vdotu(v, um2) + (vT - vp)*um2T) + umT*vp*um2p )) + + params["gamma"]*params["Ma^2"]*(ump*(ugradu(v, um2, ub) + ugradu(v, ub, um2)) + um2p*(ugradu(v, um, ub) + ugradu(v, ub, um))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradu(v, um, um2) + ugradu(v, um2, um)) - (um2T*div(v) + ugradf(v, um2T))*ump - (umT*div(v) + ugradf(v, umT))*um2p + params["1/Re"]*(umT*visc(v, um2) + um2T*visc(v, um) + visc2(v, umT, um2) + visc2(v, um2T, um)) + + vT*(params["gamma"]*params["Ma^2"]*(ump*(ugradf(um2, ubT) + ugradf(ub, um2T) + (params["gamma"] - 1.0)*(um2T*div(ub) + ubT*div(um2))) + um2p*(ugradf(um, ubT) + ugradf(ub, umT) + (params["gamma"] - 1.0)*(umT*div(ub) + ubT*div(um)))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(um, um2T) + ugradf(um2, umT) + (params["gamma"] - 1.0)*(umT*div(um2) + um2T*div(um))) - params["gamma"]*params["Ma^2"]*((params["gamma"] - 1.0)*params["1/Re"]*(umT*(visc(um2, ub) + visc(ub, um2)) + um2T*(visc(um, ub) + visc(ub, um)) + ubT*(visc(um, um2) + visc(um2, um))))) + params["gamma"]*params["1/Re"]*params["1/Pr"]*(umT*diff(vT, um2T) + um2T*diff(vT, umT) + 2.0*vT*diff(um2T, umT)) + + vp*(params["gamma"]*params["Ma^2"]*(umT*(ugradf(um2, ubp) + ugradf(ub, um2p)) + um2T*(ugradf(um, ubp) + ugradf(ub, ump)) + ubT*(ugradf(um, um2p) + ugradf(um2, ump)) + ump*(um2T*div(ub) + ubT*div(um2) - ugradf(um2, ubT) - ugradf(ub, um2T)) + um2p*(umT*div(ub) + ubT*div(um) - ugradf(um, ubT) - ugradf(ub, umT))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(umT*div(um2) + um2T*div(um) - ugradf(um, um2T) - ugradf(um2, umT))) + ) + + int2d(Th)( + params["gamma"]*params["Ma^2"]*(iomega*(dump*(vdotu(v, um) + (vT - vp)*umT) + dumT*vp*ump ) + iomega2*(ump*(vdotu(v, dum) + (vT - vp)*dumT) + umT*vp*dump )) + + params["gamma"]*params["Ma^2"]*(ump*(ugradu(v, dum, ub) + ugradu(v, ub, dum)) + dump*(ugradu(v, um, ub) + ugradu(v, ub, um))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradu(v, um, dum) + ugradu(v, dum, um)) - (dumT*div(v) + ugradf(v, dumT))*ump - (umT*div(v) + ugradf(v, umT))*dump + params["1/Re"]*(umT*visc(v, dum) + dumT*visc(v, um) + visc2(v, umT, dum) + visc2(v, dumT, um)) + + vT*(params["gamma"]*params["Ma^2"]*(ump*(ugradf(dum, ubT) + ugradf(ub, dumT) + (params["gamma"] - 1.0)*(dumT*div(ub) + ubT*div(dum))) + dump*(ugradf(um, ubT) + ugradf(ub, umT) + (params["gamma"] - 1.0)*(umT*div(ub) + ubT*div(um)))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(ugradf(um, dumT) + ugradf(dum, umT) + (params["gamma"] - 1.0)*(umT*div(dum) + dumT*div(um))) - params["gamma"]*params["Ma^2"]*((params["gamma"] - 1.0)*params["1/Re"]*(umT*(visc(dum, ub) + visc(ub, dum)) + dumT*(visc(um, ub) + visc(ub, um)) + ubT*(visc(um, dum) + visc(dum, um))))) + params["gamma"]*params["1/Re"]*params["1/Pr"]*(umT*diff(vT, dumT) + dumT*diff(vT, umT) + 2.0*vT*diff(dumT, umT)) + + vp*(params["gamma"]*params["Ma^2"]*(umT*(ugradf(dum, ubp) + ugradf(ub, dump)) + dumT*(ugradf(um, ubp) + ugradf(ub, ump)) + ubT*(ugradf(um, dump) + ugradf(dum, ump)) + ump*(dumT*div(ub) + ubT*div(dum) - ugradf(dum, ubT) - ugradf(ub, dumT)) + dump*(umT*div(ub) + ubT*div(um) - ugradf(um, ubT) - ugradf(ub, umT))) + (params["gamma"]*params["Ma^2"]*ubp + 1.0)*(umT*div(dum) + dumT*div(um) - ugradf(um, dumT) - ugradf(dum, umT))) + ) + + HomBoundaryConditions(dum); + +// TRESSIAN OPERATOR + varf vT(defu(dum), defu(v)) + = int2d(Th)( + params["gamma"]*params["Ma^2"]*(ump*(ugradu(v, um2, um3) + ugradu(v, um3, um2)) + um2p*(ugradu(v, um, um3) + ugradu(v, um3, um)) + um3p*(ugradu(v, um, um2) + ugradu(v, um2, um))) + + vT*(params["gamma"]*params["Ma^2"]*(ump*(ugradf(um2, um3T) + ugradf(um3, um2T) + (params["gamma"] - 1.0)*(um2T*div(um3) + um3T*div(um2))) + um2p*(ugradf(um, um3T) + ugradf(um3, umT) + (params["gamma"] - 1.0)*(umT*div(um3) + um3T*div(um))) + um3p*(ugradf(um, um2T) + ugradf(um2, umT) + (params["gamma"] - 1.0)*(umT*div(um2) + um2T*div(um)))) - params["gamma"]*params["Ma^2"]*((params["gamma"] - 1.0)*params["1/Re"]*(umT*(visc(um2, um3) + visc(um3, um2)) + um2T*(visc(um, um3) + visc(um3, um)) + um3T*(visc(um, um2) + visc(um2, um))))) + + vp*params["gamma"]*params["Ma^2"]*(umT*(ugradf(um2, um3p) + ugradf(um3, um2p)) + um2T*(ugradf(um, um3p) + ugradf(um3, ump)) + um3T*(ugradf(um, um2p) + ugradf(um2, ump)) + ump*(um2T*div(um3) + um3T*div(um2) - ugradf(um2, um3T) - ugradf(um3, um2T)) + um2p*(umT*div(um3) + um3T*div(um) - ugradf(um, um3T) - ugradf(um3, umT)) + um3p*(umT*div(um2) + um2T*div(um) - ugradf(um, um2T) - ugradf(um2, umT))) + ) + + int2d(Th)( + params["gamma"]*params["Ma^2"]*(ump*(ugradu(v, um2, dum) + ugradu(v, dum, um2)) + um2p*(ugradu(v, um, dum) + ugradu(v, dum, um)) + dump*(ugradu(v, um, um2) + ugradu(v, um2, um))) + + vT*(params["gamma"]*params["Ma^2"]*(ump*(ugradf(um2, dumT) + ugradf(dum, um2T) + (params["gamma"] - 1.0)*(um2T*div(dum) + dumT*div(um2))) + um2p*(ugradf(um, dumT) + ugradf(dum, umT) + (params["gamma"] - 1.0)*(umT*div(dum) + dumT*div(um))) + dump*(ugradf(um, um2T) + ugradf(um2, umT) + (params["gamma"] - 1.0)*(umT*div(um2) + um2T*div(um)))) - params["gamma"]*params["Ma^2"]*((params["gamma"] - 1.0)*params["1/Re"]*(umT*(visc(um2, dum) + visc(dum, um2)) + um2T*(visc(um, dum) + visc(dum, um)) + dumT*(visc(um, um2) + visc(um2, um))))) + + vp*params["gamma"]*params["Ma^2"]*(umT*(ugradf(um2, dump) + ugradf(dum, um2p)) + um2T*(ugradf(um, dump) + ugradf(dum, ump)) + dumT*(ugradf(um, um2p) + ugradf(um2, ump)) + ump*(um2T*div(dum) + dumT*div(um2) - ugradf(um2, dumT) - ugradf(dum, um2T)) + um2p*(umT*div(dum) + dumT*div(um) - ugradf(um, dumT) - ugradf(dum, umT)) + dump*(umT*div(um2) + um2T*div(um) - ugradf(um, um2T) - ugradf(um2, umT))) + ) + + HomBoundaryConditions(dum); + \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_structure.idp b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_structure.idp new file mode 100644 index 0000000..9a7c61c --- /dev/null +++ b/examples/douglas_jolivet_2026/eqns_douglas_jolivet_2026_structure.idp @@ -0,0 +1,132 @@ +// +// eqns_douglas_jolivet_2026_structure.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// macros +macro vdotu(v, u) (v*u + v#y*u#y + v#z*u#z) //EOM + +macro el(u) [ dx(u), + dy(u#y), + dz(u#z), + dy(u) + dx(u#y), + dz(u#y) + dy(u#z), + dx(u#z) + dz(u) + ]// EOM + +macro en(u, U) [ 0.5*(dx(u)*dx(U) + dx(u#y)*dx(U#y) + dx(u#z)*dx(U#z)), + 0.5*(dy(u)*dy(U) + dy(u#y)*dy(U#y) + dy(u#z)*dy(U#z)), + 0.5*(dz(u)*dz(U) + dz(u#y)*dz(U#y) + dz(u#z)*dz(U#z)), + dx(u)*dy(U) + dx(u#y)*dy(U#y) + dx(u#z)*dy(U#z), + dy(u)*dz(U) + dy(u#y)*dz(U#y) + dy(u#z)*dz(U#z), + dz(u)*dx(U) + dz(u#y)*dx(U#y) + dz(u#z)*dx(U#z) + ]// EOM + +macro GLE(u)(el(u) + en(u, u)) // EOM +macro dGLE(u, U) (el(u) + en(u, U) + en(U, u)) // EOM +macro ddGLE(u, U)(en(u, U) + en(U, u)) // EOM + +func D = [[1. - params["nu"], params["nu"], params["nu"], 0, 0, 0 ], + [params["nu"], 1. - params["nu"], params["nu"], 0, 0, 0 ], + [params["nu"], params["nu"], 1. - params["nu"], 0, 0, 0 ], + [0, 0, 0, 0.5 - params["nu"], 0, 0 ], + [0, 0, 0, 0, 0.5 - params["nu"], 0 ], + [0, 0, 0, 0, 0, 0.5 - params["nu"]]]; + +// Define KSP parameters + string KSPparams = "-ksp_type preonly -pc_type cholesky"; + +// Boundary conditions + macro BoundaryConditions(u, U) + on(BCxsym, u = U) + on(BCysym, u#y = U#y) + // EOM + macro HomBoundaryConditions(u) + on(((abs(int(sym(0))) % 2) == 0)*BCxsym, u = 0) + + on(((abs(int(sym(0))) % 2) != 0)*BCxsym, u#y = 0, u#z = 0) + + on(((abs(int(sym(1))) % 2) == 0)*BCysym, u#y = 0) + + on(((abs(int(sym(1))) % 2) != 0)*BCysym, u = 0, u#z = 0) + // EOM + +// Trick to deduce symmetry from the mesh and renormalize the point load +real normalizer; +{ + int[int] meshlabels = labels(Thg); + int xsym, ysym; + for (int i; i defu(um), defu(um2), defu(um3); +complex[int, int] uh(um[].n, max(1, Nh0)), umh(um[].n, max(1, 2*Nh1)); +complex eigenvalue; +if (fileext == "porb") { + ub[] = loadporb(fileroot, meshin, uh, sym0, omega, Nh0); +} +else if (fileext == "floq"){ + um[] = loadfloq(fileroot, meshin, umh, sym1, eigenvalue, sym0, omega, Nh1); + string porbfileroot, porbfilein = readbasename(workdir + filein); + parsefilename(porbfilein, porbfileroot); + ub[] = loadporb(porbfileroot, meshin, uh, sym0, omega, Nh0); +} +int isasym; +if (symmetry == "eighth" || symmetry == "Qy" || symmetry == "Qz" || symmetry == "Hyz"){ + defu(ub) = [(x >= 0)*ub(x, y, z) + (x < 0)*ub(-x, y, z), + (x >= 0)*ubY(x, y, z) + (x < 0)*ubY(-x, y, z)]; + defu(um) = [(x >= 0)*um(x, y, z) + (sym1(0) ? -1 : 1)*(x < 0)*um(-x, y, z), + (x >= 0)*umY(x, y, z) + (sym1(0) ? -1 : 1)*(x < 0)*umY(-x, y, z)]; + for (int ii = 0; ii < Nh0; ii++) { + isasym = (sym0(0)*(ii+1) % 2); + um2[] = uh(:, ii); + defu(um3) = [(x >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2(-x, y, z), + (x >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2Y(-x, y, z)]; + uh(:, ii) = um3[]; + } + for (int ii = 0; ii < Nh1; ii++) { + isasym = (sym1(0)*(ii) % 2); + um2[] = umh(:, 2*ii); + defu(um3) = [(x >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2(-x, y, z), + (x >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2Y(-x, y, z)]; + umh(:, 2*ii) = um3[]; + um2[] = umh(:, 2*ii+1); + defu(um3) = [(x >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2(-x, y, z), + (x >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(x < 0)*um2Y(-x, y, z)]; + umh(:, 2*ii+1) = um3[]; + } +} +if (symmetry == "eighth" || symmetry == "Qx" || symmetry == "Qz" || symmetry == "Hxz"){ + defu(ub) = [(y >= 0)*ub(x, y, z) + (y < 0)*ub(x, -y, z), + (y >= 0)*ubY(x, y, z) + (y < 0)*ubY(x, -y, z)]; + defu(um) = [(y >= 0)*um(x, y, z) + (sym1(1) ? -1 : 1)*(y < 0)*um(x, -y, z), + (y >= 0)*umY(x, y, z) + (sym1(1) ? -1 : 1)*(y < 0)*umY(x, -y, z)]; + for (int ii = 0; ii < Nh0; ii++) { + um2[] = uh(:, ii); + isasym = (sym0(1)*(ii+1) % 2); + defu(um3) = [(y >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2(x, -y, z), + (y >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2Y(x, -y, z)]; + uh(:, ii) = um3[]; + } + for (int ii = 0; ii < Nh1; ii++) { + isasym = (sym1(1)*(ii) % 2); + um2[] = umh(:, 2*ii); + defu(um3) = [(y >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2(x, -y, z), + (y >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2Y(x, -y, z)]; + umh(:, 2*ii) = um3[]; + um2[] = umh(:, 2*ii+1); + defu(um3) = [(y >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2(x, -y, z), + (y >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(y < 0)*um2Y(x, -y, z)]; + umh(:, 2*ii+1) = um3[]; + } +} +if (symmetry == "eighth" || symmetry == "Qx" || symmetry == "Qy" || symmetry == "Hxy"){ + defu(ub) = [(z >= 0)*ub(x, y, z) + (z < 0)*ub(x, y, -z), + (z >= 0)*ubY(x, y, z) + (z < 0)*ubY(x, y, -z)]; + defu(um) = [(z >= 0)*um(x, y, z) + (sym1(2) ? -1 : 1)*(z < 0)*um(x, y, -z), + (z >= 0)*umY(x, y, z) + (sym1(2) ? -1 : 1)*(z < 0)*umY(x, y, -z)]; + for (int ii = 0; ii < Nh0; ii++) { + isasym = (sym0(2)*(ii+1) % 2); + um2[] = uh(:, ii); + defu(um3) = [(z >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2(x, y, -z), + (z >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2Y(x, y, -z)]; + uh(:, ii) = um3[]; + } + for (int ii = 0; ii < Nh1; ii++) { + isasym = (sym1(2)*(ii) % 2); + um2[] = umh(:, 2*ii); + defu(um3) = [(z >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2(x, y, -z), + (z >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2Y(x, y, -z)]; + umh(:, 2*ii) = um3[]; + um2[] = umh(:, 2*ii+1); + defu(um3) = [(z >= 0)*um2(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2(x, y, -z), + (z >= 0)*um2Y(x, y, z) + (isasym ? -1 : 1)*(z < 0)*um2Y(x, y, -z)]; + umh(:, 2*ii+1) = um3[]; + } +} +saveporb(fileout + "start", "", meshin, sym0, omega, Nh0, true, true); +if (fileext == "floq") { +Nh0 = min(Nh0, Nh1); +ub[] += 2.0*amplitude*um[].re; +for (int ii = 0; ii < Nh0; ii++) { + uh(:, ii) += amplitude*umh(:,2*ii); + complex[int] qq = amplitude*umh(:,2*ii+1); + uh(:, ii) += conj(qq); +} +saveporb(fileout + "guess", "", meshin, sym1, omega, Nh0, true, true); +ub[] = 2.0*amplitude*um[].re; +for (int ii = 0; ii < Nh0; ii++) { + uh(:, ii) = amplitude*umh(:,2*ii); + complex[int] qq = amplitude*umh(:,2*ii+1); + uh(:, ii) += conj(qq); +} +saveporb(fileout + "branch", "", meshin, sym1, omega, Nh0, true, true); +} +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/extend_cylinder.md b/examples/douglas_jolivet_2026/extend_cylinder.md new file mode 100644 index 0000000..39ecf7e --- /dev/null +++ b/examples/douglas_jolivet_2026/extend_cylinder.md @@ -0,0 +1,59 @@ +# extend_cylinder.md +Author: Chris Douglas ([@cmdoug](https://github.com/cmdoug)) [christopher.douglas@duke.edu](mailto:christopher.douglas@duke.edu) + +MUST BE RUN WITH 1 MPI PROCESS + +```freefem +load "iovtk" +include "settings.idp" +include "macros_bifbox.idp" +// arguments +string filein = getARGV("-fi", ""); +string fileout = getARGV("-fo", ""); +real amplitude = getARGV("-amp", 1.0); + +assert(mpisize==1); + +string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension +string outfileext = parsefilename(fileout, fileout); // trim extension from output file, if given +if(fileext == "mode" || fileext == "resp" || fileext == "rslv" || fileext == "tdls" || fileext == "floq"){ + filein = readbasename(workdir + filein); + fileext = parsefilename(filein, fileroot); +} +string meshin = readmeshname(workdir + filein); // get mesh file +string meshroot, meshext = parsefilename(meshin, meshroot); +string symmetry = meshroot(meshroot.rfind("_")+1:meshroot.length-1); // get file extension + +// Load mesh, make FE basis +Th = readmeshN(workdir + meshin); +meshin = meshroot + "full." + meshext; // add mesh extension for full domain +mesh Th2 = movemesh(Th,[x,-y]); +Th = Th + Th2; +Th = change(Th, rmInternalEdges=1); +savemesh(Th, workdir + meshin); +Thg = Th; +DmeshCreate(Th); +restu = restrict(XMh, XMhg, n2o); +XMh defu(ub); +XMh defu(um), defu(uma); + +int Nh = getARGV("-Nh", 1); +complex[int, int] uh(ub[].n, Nh); +real omega; +ub[] = loadporb(fileroot, meshin, uh, sym, omega, Nh); + +defu(ub) = [(y >= 0)*ub(x, y) + (y < 0)*ub(x, -y), + (y >= 0)*uby(x, y) - (y < 0)*uby(x, -y), + (y >= 0)*ubT(x, y) + (y < 0)*ubT(x, -y), + (y >= 0)*ubp(x, y) + (y < 0)*ubp(x, -y)]; +for (int ii = 0; ii < Nh; ii++) { + uma[] = 2.*uh(:, ii).re; + int isasym = sym(0)*(ii+1); + defu(um) = [(y >= 0)*uma(x, y) + (isasym ? -1 : 1)*(y < 0)*uma(x, -y), + (y >= 0)*umay(x, y) - (isasym ? -1 : 1)*(y < 0)*umay(x, -y), + (y >= 0)*umaT(x, y) + (isasym ? -1 : 1)*(y < 0)*umaT(x, -y), + (y >= 0)*umap(x, y) + (isasym ? -1 : 1)*(y < 0)*umap(x, -y)]; + ub[] += um[]; +} +savebase(fileout, "", meshin, true, true); +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/extend_structure.md b/examples/douglas_jolivet_2026/extend_structure.md new file mode 100644 index 0000000..d371908 --- /dev/null +++ b/examples/douglas_jolivet_2026/extend_structure.md @@ -0,0 +1,76 @@ +# extend_structure.md +Author: Chris Douglas ([@cmdoug](https://github.com/cmdoug)) [christopher.douglas@duke.edu](mailto:christopher.douglas@duke.edu) + +MUST BE RUN WITH 1 MPI PROCESS + +```freefem +load "iovtk" +include "settings.idp" +include "macros_bifbox.idp" +// arguments +string filein = getARGV("-fi", ""); +string fileout = getARGV("-fo", ""); +real amplitude = getARGV("-amp", 1.0); + +assert(mpisize==1); + +string fileroot, fileext = parsefilename(filein, fileroot); //extract file name and extension +string outfileext = parsefilename(fileout, fileout); // trim extension from output file, if given +if(fileext == "mode" || fileext == "resp" || fileext == "rslv" || fileext == "tdls" || fileext == "floq"){ + filein = readbasename(workdir + filein); + fileext = parsefilename(filein, fileroot); +} +string meshin = readmeshname(workdir + filein); // get mesh file +string meshroot, meshext = parsefilename(meshin, meshroot); +string symmetry = meshroot(meshroot.rfind("_")+1:meshroot.length-1); // get file extension +meshin = meshroot(0:meshroot.rfind("_")-1); // get file root less symmetry and file extension +meshin = meshin + "_A." + meshext; // add mesh extension for full domain + +// Load mesh, make FE basis +Th = readmeshN(workdir + meshin); +Thg = Th; +DmeshCreate(Th); +restu = restrict(XMh, XMhg, n2o); +XMh defu(ub); +XMh defu(um), defu(uma); + +if (fileext == "base") ub[] = loadbase(fileroot, meshin); +else if(fileext == "fold") { + real[string] alpha; + real beta; + ub[] = loadfold(fileroot, meshin, um[], uma[], alpha, beta); +} +else if (fileext == "hopf") { + complex eigenvalue; + real omega; + complex[string] alpha; + complex beta; + complex[int] qm(um[].n), qma(um[].n); + ub[] = loadhopf(fileroot, meshin, qm, qma, sym, omega, alpha, beta); + um[] = 2.*qm.re; + uma[] = 2.*qma.re; +} +if (symmetry == "S" || symmetry == "Hx"){ + defu(ub) = [(x <= 0)*ub(x, y, z) - (x > 0)*ub(-x, y, z), + (x <= 0)*uby(x, y, z) + (x > 0)*uby(-x, y, z), + (x <= 0)*ubz(x, y, z) + (x > 0)*ubz(-x, y, z)]; + defu(um) = [(x <= 0)*um(x, y, z) - (sym(0) ? -1 : 1)*(x > 0)*um(-x, y, z), + (x <= 0)*umy(x, y, z) + (sym(0) ? -1 : 1)*(x > 0)*umy(-x, y, z), + (x <= 0)*umz(x, y, z) + (sym(0) ? -1 : 1)*(x > 0)*umz(-x, y, z)]; +} + +if (symmetry == "S" || symmetry == "Hy"){ + defu(ub) = [(y <= 0)*ub(x, y, z) + (y > 0)*ub(x, -y, z), + (y <= 0)*uby(x, y, z) - (y > 0)*uby(x, -y, z), + (y <= 0)*ubz(x, y, z) + (y > 0)*ubz(x, -y, z)]; + defu(um) = [(y <= 0)*um(x, y, z) + (sym(1) ? -1 : 1)*(y > 0)*um(x, -y, z), + (y <= 0)*umy(x, y, z) - (sym(1) ? -1 : 1)*(y > 0)*umy(x, -y, z), + (y <= 0)*umz(x, y, z) + (sym(1) ? -1 : 1)*(y > 0)*umz(x, -y, z)]; +} + +savebase(fileout + "start", "", meshin, true, true); +ub[] += amplitude*um[]; +savebase(fileout + "guess", "", meshin, true, true); +ub[] = amplitude*um[]; +savebase(fileout + "branch", "", meshin, true, true); +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/mesh_brusselator.md b/examples/douglas_jolivet_2026/mesh_brusselator.md new file mode 100644 index 0000000..6879882 --- /dev/null +++ b/examples/douglas_jolivet_2026/mesh_brusselator.md @@ -0,0 +1,86 @@ +# mesh_brusselator.md +Author: Chris Douglas ([@cmdoug](https://github.com/cmdoug)) [christopher.douglas@duke.edu](mailto:christopher.douglas@duke.edu) + +This file is used to create initial meshes for the example of the 3-D Brusselator from [C. Douglas, P. Jolivet. (2026)]. + +```freefem +assert(mpisize == 1); // Must be run with 1 processor +include "settings.idp" // so that boundary labels match the convention in the settings file +int n2 = ceil(getARGV("-n", 10)/2); +string meshout = getARGV("-mo", "cube"); // mesh filename +if(meshout.rfind(".mesh") > 0) meshout = meshout(0:meshout.length - 6); // add extension if not provided + +// ASSEMBLE CUBE MESH USING REFLECTIONS +// base 1/8 mesh +mesh3 Th8 = cube(n2, n2, n2, flags = 5); // create mesh +Th8 = movemesh(Th8, [x/2, y/2, z/2]); +int[int] labs = [1, BCsY, 2, BCnX, 3, BCnY, 4, BCsX, 5, BCsZ, 6, BCnZ]; +Th8 = change(Th8, label = labs); +int[int] meshlabels = labels(Th8); +cout << "\tMesh: " << Th8.nv << " vertices, " << Th8.nt << " elements, " << Th8.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_eighth.mesh'." << endl; +savemesh(Th8, meshout + "_eighth.mesh"); +// 1/4 mesh -- reflected across x +mesh3 Th4x = movemesh(Th8, [-x, y, z], orientation = -1); // move mesh to center +Th4x = Th8 + Th4x; +labs = [1, BCsY, 2, BCnX, 3, BCnY, 4, BCdX, 5, BCsZ, 6, BCnZ]; +Th4x = change(Th4x, label = labs, rmInternalFaces = true); +meshlabels = labels(Th4x); +cout << "\tMesh: " << Th4x.nv << " vertices, " << Th4x.nt << " elements, " << Th4x.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Qx.mesh'." << endl; +savemesh(Th4x, meshout + "_Qx.mesh"); +// 1/4 mesh -- reflected across y +mesh3 Th4y = movemesh(Th8, [x, -y, z], orientation = -1); // move mesh to center +Th4y = Th8 + Th4y; +labs = [1, BCdY, 2, BCnX, 3, BCnY, 4, BCsX, 5, BCsZ, 6, BCnZ]; +Th4y = change(Th4y, label = labs, rmInternalFaces = true); +meshlabels = labels(Th4y); +cout << "\tMesh: " << Th4y.nv << " vertices, " << Th4y.nt << " elements, " << Th4y.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Qy.mesh'." << endl; +savemesh(Th4y, meshout + "_Qy.mesh"); +// 1/4 mesh -- reflected across z +mesh3 Th4z = movemesh(Th8, [x, y, -z], orientation = -1); // move mesh to center +Th4z = Th8 + Th4z; +labs = [1, BCsY, 2, BCnX, 3, BCnY, 4, BCsX, 5, BCdZ, 6, BCnZ]; +Th4z = change(Th4z, label = labs, rmInternalFaces = true); +meshlabels = labels(Th4z); +cout << "\tMesh: " << Th4z.nv << " vertices, " << Th4z.nt << " elements, " << Th4z.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Qz.mesh'." << endl; +savemesh(Th4z, meshout + "_Qz.mesh"); +// 1/2 mesh -- reflected across x,y +mesh3 Th2xy = movemesh(Th4x, [x, -y, z], orientation = -1); // move mesh to center +Th2xy = Th4x + Th2xy; +labs = [1, BCdY, 2, BCnX, 3, BCnY, 4, BCdX, 5, BCsZ, 6, BCnZ]; +Th2xy = change(Th2xy, label = labs, rmInternalFaces = true); +meshlabels = labels(Th2xy); +cout << "\tMesh: " << Th2xy.nv << " vertices, " << Th2xy.nt << " elements, " << Th2xy.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Hxy.mesh'." << endl; +savemesh(Th2xy, meshout + "_Hxy.mesh"); +// 1/2 mesh -- reflected across x,z +mesh3 Th2xz = movemesh(Th4x, [x, y, -z], orientation = -1); // move mesh to center +Th2xz = Th4x + Th2xz; +labs = [1, BCsY, 2, BCnX, 3, BCnY, 4, BCdX, 5, BCdZ, 6, BCnZ]; +Th2xz = change(Th2xz, label = labs, rmInternalFaces = true); +meshlabels = labels(Th2xz); +cout << "\tMesh: " << Th2xz.nv << " vertices, " << Th2xz.nt << " elements, " << Th2xz.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Hxz.mesh'." << endl; +savemesh(Th2xz, meshout + "_Hxz.mesh"); +// 1/2 mesh -- reflected across y,z +mesh3 Th2yz = movemesh(Th4y, [x, y, -z], orientation = -1); // move mesh to center +Th2yz = Th4y + Th2yz; +labs = [1, BCdY, 2, BCnX, 3, BCnY, 4, BCsX, 5, BCdZ, 6, BCnZ]; +Th2yz = change(Th2yz, label = labs, rmInternalFaces = true); +meshlabels = labels(Th2yz); +cout << "\tMesh: " << Th2yz.nv << " vertices, " << Th2yz.nt << " elements, " << Th2yz.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_Hyz.mesh'." << endl; +savemesh(Th2yz, meshout + "_Hyz.mesh"); +// full mesh -- reflected across x,y,z +mesh3 Thg = movemesh(Th2xy, [x, y, -z], orientation = -1); // move mesh to center +Thg = Thg + Th2xy; +labs = [1, BCdY, 2, BCnX, 3, BCnY, 4, BCdX, 5, BCdZ, 6, BCnZ]; +Thg = change(Thg, label = labs, rmInternalFaces = true); +meshlabels = labels(Thg); +cout << "\tMesh: " << Thg.nv << " vertices, " << Thg.nt << " elements, " << Thg.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "_full.mesh'." << endl; +savemesh(Thg, meshout + "_full.mesh"); +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/mesh_cylinder.md b/examples/douglas_jolivet_2026/mesh_cylinder.md new file mode 100644 index 0000000..f0d9053 --- /dev/null +++ b/examples/douglas_jolivet_2026/mesh_cylinder.md @@ -0,0 +1,37 @@ +# mesh_cylinder.md +Author: Chris Douglas ([@cmdoug](https://github.com/cmdoug)) [christopher.douglas@duke.edu](mailto:christopher.douglas@duke.edu) + +This file is used to create an initial mesh for the example of 2-D flow past a cylinder from [C. Douglas, P. Jolivet. (2026)]. + +```freefem +assert(mpisize == 1); // Must be run with 1 processor +include "settings.idp" // so that boundary labels match the convention in the settings file +real n0 = 50; +real n1 = 2.0; +real n2 = 0.25; +real xmin = -100.0; +real xmax = 150.0; +real rmax = 100.0; +real rcyl = 0.5; +string meshout = getARGV("-mo", "cylinder.msh"); // mesh filename +if(meshout.rfind(".msh") < 0) meshout = meshout + ".msh"; // add extension if not provided + +// Define borders +// o------------------6-------------------o +// | | +// 1 5 +// | ╭--3--╮ | +// o-----2-----o 0 o---------4----------o +border C01(t=0, 1){x=xmin; y=rmax*(1-t); label=BCin;} +border C02(t=0, 1){x=xmin+(abs(xmin)-rcyl)*t; y=0; label=BCsym;} +border C03(t=0, 1){x=-rcyl*cos(pi*t); y=rcyl*sin(pi*t); label=BCwall;} +border C04(t=0, 1){x=rcyl+(xmax-rcyl)*t; y=0; label=BCsym;} +border C05(t=0, 1){x=xmax; y=rmax*t; label=BCout;} +border C06(t=0, 1){x=xmax+(xmin-xmax)*t; y=rmax; label=BClat;} +// Assemble mesh +mesh Thg = buildmesh(C01(n2*rmax) + C02(n1*(abs(xmin)-rcyl)) + C03(n0*pi*rcyl) + C04(n1*(xmax-rcyl)) + C05(n2*rmax) + C06(n2*(xmax-xmin))); +int[int] meshlabels = labels(Thg); +cout << "\tMesh: " << Thg.nv << " vertices, " << Thg.nt << " elements, " << Thg.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; +cout << " Saving mesh '" + meshout + "'." << endl; +savemesh(Thg, meshout); +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/mesh_structure.md b/examples/douglas_jolivet_2026/mesh_structure.md new file mode 100644 index 0000000..6f22ab4 --- /dev/null +++ b/examples/douglas_jolivet_2026/mesh_structure.md @@ -0,0 +1,67 @@ +# mesh_structure.md +Author: Chris Douglas ([@cmdoug](https://github.com/cmdoug)) [christopher.douglas@duke.edu](mailto:christopher.douglas@duke.edu) + +This file is used to create initial meshes for the example of the buckling 3-D structure from [C. Douglas, P. Jolivet. (2026)]. + +```freefem +assert(mpisize == 1); // Must be run with 1 processor +include "settings.idp" // so that boundary labels match the convention in the settings file +real R = 2.540; //[m] +real t = 0.00635; //[m] +real beta = 0.1; //[rad] +int n0 = 2*ceil(getARGV("-n0", 40)/2); //must be even for discrete symmetry and point loading +int n1 = 2*ceil(getARGV("-n1", 2)/2); //must be even for pinned BC +int flag = getARGV("-flag", 1); +string meshroot, meshout = getARGV("-mo", "shell.mesh"); // mesh filename + +func string parsefilename(string & filename, string & fileroot){ + string fileext; + if(filename.rfind(".") > 0){ // filename includes extension + fileext = filename(filename.rfind(".")+1:filename.length-1); // get file extension + fileroot = filename(0:filename.rfind(".")-1); // get file root + } + else { + fileroot = filename; + fileext = ""; + } + return fileext; +} + +if(meshout.rfind(".mesh") < 0) meshout = meshout + ".mesh"; // add extension if not provided +string meshext = parsefilename(meshout, meshroot); + +string[int] symlist = ["S", "Hx", "Hy", "Rz", "A"]; + +int[int] rup = [0, BCtop], rdown = [0, BCfree], rmid(8); +mesh Th = square(n0/2, n0/2, [beta*(x - 1.), beta*(y - 1.)], flags = flag, region = 0); +int[int] regnum = [0, 0]; +mesh3 Th3 = buildlayers(Th, n1, zbound = [-t/2., t/2.], region = regnum, labelup = rup, labeldown = rdown); +for (int ii = 0; ii < 5; ii++) { + string meshsym = symlist[ii]; + if (meshsym == "Rz") regnum = [0, BCpsym]; + else regnum = [0, 0]; + mesh3 Thg = movemesh(Th3, [R*x, (R+z)*sin(y), (R+z)*cos(y) - R*cos(beta)], region = regnum); + if (meshsym == "S") rmid = [1, BCpin, 2, BCxsym, 3, BCysym, 4, BCfree]; + else if (meshsym == "Hx") { + mesh3 Thy = movemesh(Thg, [x, -y, z], orientation = -1); + Thg = Thg + Thy; + rmid = [1, BCpin, 2, BCxsym, 3, BCfree, 4, BCfree]; + } + else if (meshsym == "Hy" || meshsym == "A" || meshsym == "Rz") { + mesh3 Thx = movemesh(Thg, [-x, y, z], orientation = -1); + Thg = Thg + Thx; + if (meshsym == "A" || meshsym == "Rz"){ + Thx = movemesh(Thg, [x, -y, z], orientation = -1); + Thg = Thg + Thx; + rmid = [1, BCpin, 2, BCfree, 3, BCfree, 4, BCfree]; + } + else rmid = [1, BCpin, 2, BCfree, 3, BCysym, 4, BCfree]; + } + else exit(1); + Thg = change(Thg, label = rmid, rmInternalFaces = true); + int[int] meshlabels = labels(Thg); + cout << "\tMesh: " << Thg.nv << " vertices, " << Thg.nt << " elements, " << Thg.nbe << " boundary elements, " << meshlabels.n << " labeled boundaries." << endl; + cout << " Saving mesh '" + meshroot + "_" + meshsym + "." + meshext + "'." << endl; + savemesh(Thg, meshroot + "_" + meshsym + "." + meshext); +} +``` \ No newline at end of file diff --git a/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_brusselator.idp b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_brusselator.idp new file mode 100644 index 0000000..967fc2a --- /dev/null +++ b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_brusselator.idp @@ -0,0 +1,59 @@ +// +// settings_douglas_jolivet_2026_brusselator.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// Set dimension macro for 'macro_ddm.idp' + macro dimension()3//EOM + macro meshtype()V//EOM +// Declare that a cubic nonlinear term exists and its "Tressian" operator is defined + macro cubic()//EOM +// Load hpddm macros + include "macro_ddm.idp" + verbosity = getARGV("-v",0); +// Define parameter and monitor names + string[int] paramnames = ["A", "B", "Dx", "Dy", "1/L^2"]; // set parameter names + string[int] monitornames = ["L", "Xavg", "Yavg", "L2diff"]; // set monitor names +// Declare symmetries + real[int] sym(3); +// Define state vector and FE space + macro defu(u)[u, u#Y]//EOM + macro initu(i)[i, i]//EOM + func Pk = [P2, P2]; +// Define forcing vector and FE space (for resolvent analysis) + macro deff(f)[f, f#Y]//EOM + macro initf(i)[i, i]//EOM + func Pkf = [P2, P2]; +// Define quantities for mesh adaptation and plotting in Paraview + macro adaptu(u)u, u#Y//EOM + macro adaptf(f)f, f#Y//EOM +// Name and order for real Paraview outputs + string ParaviewDataName = "X Y"; + string ParaviewDataNamef = ParaviewDataName; + int[int] ParaviewOrder = [1, 1]; + int[int] ParaviewOrderf = ParaviewOrder; +// Name and order for complex Paraview outputs + string ParaviewDataNamec = "X_r Y_r X_i Y_i"; + string ParaviewDataNamefc = ParaviewDataNamec; + int[int] ParaviewOrderc = [ParaviewOrder, ParaviewOrder]; + int[int] ParaviewOrderfc = ParaviewOrderc; + // Initial conditions (if no file) + macro InitialConditions()[params["A"], params["B"]/params["A"]]//EOM +// Boundary labels + int BCnull = 0; + int BCdY = 1; + int BCnX = 2; + int BCnY = 3; + int BCdX = 4; + int BCdZ = 5; + int BCnZ = 6; + int BCsX = 7; + int BCsY = 8; + int BCsZ = 9; +// Define solution monitors to extract: + macro getmonitors(){ + monitors["L"] = 1.0/sqrt(params["1/L^2"]); + monitors["Xavg"] = real(int3d(Thg)(ubg)/int3d(Thg)(1.0)); + monitors["Yavg"] = real(int3d(Thg)(ubgY)/int3d(Thg)(1.0)); + monitors["L2diff"] = real(int3d(Thg)(sqrt((ubg - params["A"])^2.0 + (ubgY - params["B"]/params["A"])^2.0))/int3d(Thg)(1.0)); + }// EOM diff --git a/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_cylinder.idp b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_cylinder.idp new file mode 100644 index 0000000..ecf0658 --- /dev/null +++ b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_cylinder.idp @@ -0,0 +1,60 @@ +// +// settings_douglas_jolivet_2026_cylinder.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// Set dimension macro for 'macro_ddm.idp' (2) + macro dimension()2//EOM +// Declare that a cubic nonlinear term exists and its "Tressian" operator is defined + macro cubic()//EOM +// Load hpddm macros + include "macro_ddm.idp" + verbosity = getARGV("-v",0); +// Define parameter and monitor names + string[int] paramnames = ["1/Re", "1/Pr", "Ma^2", "gamma"]; // set parameter names + string[int] monitornames = ["uxmin", "Fx", "Fy", "Re", "Ma"]; // set monitor names +// Declare symmetries + real[int] sym(1); +// Define state vector and FE space + macro defu(u)[u, u#y, u#T, u#p]//EOM + macro initu(i)[i, i, i, i]//EOM + func Pk = [P2, P2, P2, P1]; +// Define forcing vector and FE space (for resolvent analysis) + macro deff(f)[f, f#y, f#T, f#p]//EOM + macro initf(i)[i, i, i, i]//EOM + func Pkf = [P2, P2, P2, P1]; +// Define quantities for mesh adaptation and plotting in Paraview + macro adaptu(u)[u, u#y], u#T, u#p//EOM + macro adaptf(f)[f, f#y], u#T, u#p//EOM + macro paraviewu(u)[u, u#y, 0*u], u#T, u#p//EOM + macro paraviewf(f)[f, f#y, 0*f], f#T, f#p//EOM +// Name and order for real Paraview outputs + string ParaviewDataName = "velocity temperature pressure"; + string ParaviewDataNamef = "momentum energy mass"; + int[int] ParaviewOrder = [1, 1, 1]; + int[int] ParaviewOrderf = [1, 1, 1]; +// Name and order for complex Paraview outputs + string ParaviewDataNamec = "velocity_r temperature_r pressure_r velocity_i temperature_i pressure_i"; + string ParaviewDataNamefc = "momentum_r energy_r mass_r momentum_i energy_i mass_i"; + int[int] ParaviewOrderc = [ParaviewOrder, ParaviewOrder]; + int[int] ParaviewOrderfc = [ParaviewOrderf, ParaviewOrderf]; + // Initial conditions (if no file) + macro InitialConditions()[1, 0, 1, 0]//EOM +// Boundary labels + int BCnull = 0; + int BCsym = 1; + int BCout = 2; + int BCin = 3; + int BCwall = 4; + int BClat = 5; +// Define solution monitors to extract: here x and y forces on BCwall + macro getmonitors(){ + fespace Xh1(Thg, P2); + Xh1 UX = ubg; + monitors["uxmin"] = UX[].min; + monitors["Fx"] = int1d(Thg, BCwall)(ubgp*N.x - params["1/Re"]*(N.x*(dx(ubg) - dy(ubgy)) + N.y*(dx(ubgy) + dy(ubg)))); + monitors["Fy"] = int1d(Thg, BCwall)(ubgp*N.y - params["1/Re"]*(N.x*(dx(ubgy) + dy(ubg)) + N.y*(dy(ubgy) - dx(ubg)))); + monitors["Re"] = 1.0/params["1/Re"]; + monitors["Ma"] = sqrt(params["Ma^2"]); + } + // EOM diff --git a/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp new file mode 100644 index 0000000..fa9d799 --- /dev/null +++ b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp @@ -0,0 +1,60 @@ +// +// settings_douglas_jolivet_2026_structure.idp +// Chris Douglas +// christopher.douglas@duke.edu +// +// Set dimension macro for 'macro_ddm.idp' + macro dimension()3//EOM + macro meshtype()V//EOM +// Declare that a cubic nonlinear term exists and its "Tressian" operator is defined + macro cubic()//EOM +// Load hpddm macros + include "macro_ddm.idp" + verbosity = getARGV("-v", 0); +// Define parameter and monitor names +// NOTE: The load parameter "P" has been reduced here to eliminate Young's modulus, E. +// It may be related to the "actual" load as follows: P_actual = P*E/(1+nu)/(1-2*nu) + string[int] paramnames = ["P", "nu"]; // set parameter names + string[int] monitornames = ["dispX", "dispY", "dispZ"]; // set monitor names +// Declare symmetries + real[int] sym(2); +// Define state vector and FE space + macro defu(u)[u, u#y, u#z]//EOM + macro initu(i)[i, i, i]//EOM + func Pk = [P2, P2, P2]; +// Boundary labels + int BCnull = 0; + int BCpin = 1; + int BCfree = 5; + int BCtop = 6; + int BCxsym = 7; + int BCysym = 8; + int BCpsym = 9; +// Define forcing vector and FE space (for resolvent analysis) + macro deff(f)[f, f#y, f#z]//EOM + macro initf(i)[i, i, i]//EOM + func Pkf = [P2, P2, P2]; +// Define quantities for mesh adaptation and plotting in Paraview + macro adaptu(u)[u, u#y, u#z]//EOM + macro adaptf(f)[f, f#y, f#z]//EOM +// Name and order for real Paraview outputs + string ParaviewDataName = "disp"; + string ParaviewDataNamef = ParaviewDataName; + int[int] ParaviewOrder = [1]; + int[int] ParaviewOrderf = ParaviewOrder; +// Name and order for complex Paraview outputs + string ParaviewDataNamec = ParaviewDataName + "_r " + ParaviewDataName + "_i"; + string ParaviewDataNamefc = ParaviewDataNamec; + int[int] ParaviewOrderc = [ParaviewOrder, ParaviewOrder]; + int[int] ParaviewOrderfc = ParaviewOrderc; + // Initial conditions (if no file) + macro InitialConditions()initu(0)//EOM +// coordinate mapping macros --------------------------------------------------- + macro coordinatetransform(U) x + real(U), y + real(U#y), z + real(U#z) // EOM +// Define solution monitors to extract: + macro getmonitors(){ + real z0 = 2.54*(1. - cos(0.1)); + monitors["dispX"] = real(1000.*ubg(0.,0.,z0)); + monitors["dispY"] = real(1000.*ubgy(0.,0.,z0)); + monitors["dispZ"] = -real(1000.*ubgz(0.,0.,z0)); + }// EOM diff --git a/examples/douglas_jolivet_2026/structure/.gitkeep b/examples/douglas_jolivet_2026/structure/.gitkeep new file mode 100644 index 0000000..e69de29 From 667e62fe877d063cde8655aab6f54e0319a134d3 Mon Sep 17 00:00:00 2001 From: cmdoug <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 09:09:39 -0400 Subject: [PATCH 2/7] Adds citation information --- examples/douglas_jolivet_2026/README_brusselator.md | 8 ++++++-- examples/douglas_jolivet_2026/README_cylinder.md | 10 +++++++--- examples/douglas_jolivet_2026/README_structure.md | 8 ++++++-- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/examples/douglas_jolivet_2026/README_brusselator.md b/examples/douglas_jolivet_2026/README_brusselator.md index 6a8089e..b3d705b 100644 --- a/examples/douglas_jolivet_2026/README_brusselator.md +++ b/examples/douglas_jolivet_2026/README_brusselator.md @@ -1,10 +1,14 @@ -# Examples: Douglas & Jolivet (202x) +# Examples: Douglas & Jolivet (2026) This file shows an example `ff-bifbox` workflow for reproducing the results of the study: ```tex @article{douglas_jolivet_2026, title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, - author={Douglas, Christopher M. and Jolivet, Pierre R.}, + author={Douglas, Christopher M. and Jolivet, Pierre}, year={2026}, + journal={Computer Physics Communications}, + publisher={Elsevier}, + notes={Accepted manuscript}, + doi={10.48550/arXiv.2509.18429} } ``` The commands below illustrate how to perform a bifurcation analysis of the Brusselator in 3-D using `ff-bifbox`. diff --git a/examples/douglas_jolivet_2026/README_cylinder.md b/examples/douglas_jolivet_2026/README_cylinder.md index e154f14..14b7f46 100644 --- a/examples/douglas_jolivet_2026/README_cylinder.md +++ b/examples/douglas_jolivet_2026/README_cylinder.md @@ -1,10 +1,14 @@ -# Compressible Flow Example: -This file shows an example `ff-bifbox` workflow for reproducing the results of the paper: +# Examples: Douglas & Jolivet (2026) +This file shows an example `ff-bifbox` workflow for reproducing the results of the study: ```tex @article{douglas_jolivet_2026, title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, - author={Douglas, Christopher M. and Jolivet, Pierre R.}, + author={Douglas, Christopher M. and Jolivet, Pierre}, year={2026}, + journal={Computer Physics Communications}, + publisher={Elsevier}, + notes={Accepted manuscript}, + doi={10.48550/arXiv.2509.18429} } ``` The commands below illustrate how to analyze a 2-D compressible flow past a cylinder using `ff-bifbox`. diff --git a/examples/douglas_jolivet_2026/README_structure.md b/examples/douglas_jolivet_2026/README_structure.md index 8fb3ad9..f28730d 100644 --- a/examples/douglas_jolivet_2026/README_structure.md +++ b/examples/douglas_jolivet_2026/README_structure.md @@ -1,10 +1,14 @@ -# Examples: Douglas & Jolivet (202x) +# Examples: Douglas & Jolivet (2026) This file shows an example `ff-bifbox` workflow for reproducing the results of the study: ```tex @article{douglas_jolivet_2026, title={ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs}, - author={Douglas, Christopher M. and Jolivet, Pierre R.}, + author={Douglas, Christopher M. and Jolivet, Pierre}, year={2026}, + journal={Computer Physics Communications}, + publisher={Elsevier}, + notes={Accepted manuscript}, + doi={10.48550/arXiv.2509.18429} } ``` The commands below illustrate how to perform a bifurcation analysis of buckling 3-D structure using `ff-bifbox`. From 39e35ba9ef3c4b9ca0781f67231319c8ec1be7ef Mon Sep 17 00:00:00 2001 From: cmdoug <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 09:13:26 -0400 Subject: [PATCH 3/7] add journal to example readme --- examples/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/README.md b/examples/README.md index 5041e8c..aa414d0 100644 --- a/examples/README.md +++ b/examples/README.md @@ -27,7 +27,7 @@ This folder contains example implementations that reproduce selected results fro | [douglas_etal_2021](./douglas_etal_2021) | bifurcation analysis of an incompressible axisymmetric swirling jet with axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2021). Nonlinear dynamics of fully developed swirling jets. _Journal of Fluid Mechanics_, 924, A14. doi:[10.1017/jfm.2021.615](https://doi.org/10.1017/jfm.2021.615). | | [douglas_etal_2022](./douglas_etal_2022) | bifurcation analysis of an incompressible axisymmetric swirling annular jet with coordinate transformation and axisymmetry breaking | Douglas, C. M., Emerson, B. L., & Lieuwen, T. C. (2022). Dynamics and bifurcations of laminar annular swirling and non-swirling jets. _Journal of Fluid Mechanics_, 943, A35. doi:[10.1017/jfm.2022.453](https://doi.org/10.1017/jfm.2022.453). | | [douglas_etal_2023](./douglas_etal_2023) | bifurcation analysis of a reacting weakly-compressible axisymmetric jet with axisymmetry breaking | Douglas, C. M., Polifke, W., & Lesshafft, L. (2023). Flash-back, blow-off, and symmetry breaking of premixed conical flames. _Combustion and Flame_, 258(2), 113060. doi:[10.1016/j.combustflame.2023.113060](https://doi.org/10.1016/j.combustflame.2023.113060). | -| [douglas_jolivet_2026](./douglas_jolivet_2026) | bifurcation analysis of a 3-D Brusselator system, a 3-D plate buckling system, and a 2-D compressible Navier--Stokes system | Douglas, C. M., & Jolivet, P. (2026). ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs. doi:[10.48550/arXiv.2509.18429](https://doi.org/10.48550/arXiv.2509.18429). | +| [douglas_jolivet_2026](./douglas_jolivet_2026) | bifurcation analysis of a 3-D Brusselator system, a 3-D plate buckling system, and a 2-D compressible Navier--Stokes system | Douglas, C. M., & Jolivet, P. (2026). ff-bifbox: A scalable, open-source toolbox for bifurcation analysis of nonlinear PDEs. _Computer Physics Communications_, (accepted manuscript). doi:[10.48550/arXiv.2509.18429](https://doi.org/10.48550/arXiv.2509.18429). | | [fani_etal_2018](./fani_etal_2018) | bifurcation analysis of a compressible 2-D flow over a cylinder with reflective symmetry breaking | Fani, A., Citro, V., Giannetti, F., & Auteri, F. (2018). Computation of the bluff-body sound generation by a self-consistent mean flow formulation. _Physics of Fluids_, 30(3), 036102. doi:[10.1063/1.4997536](https://doi.org/10.1063/1.4997536). | | [FK_problem](./FK_problem) | bifurcation analysis of Frank–Kamenetskii thermal runaway (explosion) in a cylinderical vessel | \Zeldovich, I. B., Barenblatt, G. I., Librovich, V. B., & Makhviladze, G. M. _Mathematical theory of combustion and explosions_, 1985. | | [garnaud_2012](./garnaud_2012) | resolvent analysis of an incompressible axisymmetric laminar jet | Garnaud, X. (2012). _Modes, transient dynamics and forced response of circular jets_. [Doctoral dissertation, Ecole Polytechnique], HAL ID:[tel-00740133](https://theses.hal.science/tel-00740133). | From 6b8aaa0da0a79334f61e7cb5287ebf32b2a84415 Mon Sep 17 00:00:00 2001 From: Chris Douglas <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 12:50:52 -0400 Subject: [PATCH 4/7] whitespace Co-authored-by: Pierre Jolivet --- examples/douglas_jolivet_2026/README_cylinder.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/douglas_jolivet_2026/README_cylinder.md b/examples/douglas_jolivet_2026/README_cylinder.md index 14b7f46..1a88b5c 100644 --- a/examples/douglas_jolivet_2026/README_cylinder.md +++ b/examples/douglas_jolivet_2026/README_cylinder.md @@ -86,7 +86,7 @@ ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi cylinder_21.hopf -fo 7. Continue the branches of periodic solutions emanating from the Hopf points along $Re$ and $Ma$. ```sh for Ma in 0 4 6 8; do -ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma".hopf -fo cylinder_Ma0p"$Ma" -mo cyl_Ma0p"$Ma" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param 1/Re -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 + ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma".hopf -fo cylinder_Ma0p"$Ma" -mo cyl_Ma0p"$Ma" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param 1/Re -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 done for Re in 50 70 90; do ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Re"$Re".hopf -fo cylinder_Re"$Re" -mo cyl_Re"$Re" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param Ma^2 -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 From bd2f9e04df826a34c38df5ad16f10177fb0cc92c Mon Sep 17 00:00:00 2001 From: Chris Douglas <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 12:51:03 -0400 Subject: [PATCH 5/7] whitespace Co-authored-by: Pierre Jolivet --- examples/douglas_jolivet_2026/README_cylinder.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/douglas_jolivet_2026/README_cylinder.md b/examples/douglas_jolivet_2026/README_cylinder.md index 1a88b5c..2d7d0b9 100644 --- a/examples/douglas_jolivet_2026/README_cylinder.md +++ b/examples/douglas_jolivet_2026/README_cylinder.md @@ -89,6 +89,6 @@ for Ma in 0 4 6 8; do ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Ma0p"$Ma".hopf -fo cylinder_Ma0p"$Ma" -mo cyl_Ma0p"$Ma" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param 1/Re -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 done for Re in 50 70 90; do -ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Re"$Re".hopf -fo cylinder_Re"$Re" -mo cyl_Re"$Re" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param Ma^2 -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 + ff-mpirun -np $nproc porbcontinue.md -v 0 -dir $workdir -fi cylinder_Re"$Re".hopf -fo cylinder_Re"$Re" -mo cyl_Re"$Re" -thetamax 0.1 -hmin 1e-5 -hmax 2 -param Ma^2 -h0 1 -scount 4 -maxcount -1 -paramtarget 0.01 -Nh 3 -fieldsplit_0_fieldsplit_0_mat_mumps_icntl_35 1 -fieldsplit_0_fieldsplit_0_mat_mumps_cntl_7 1.0e-8 done ``` \ No newline at end of file From facce932d93941dc2f839e5a3fe0c394e70a116d Mon Sep 17 00:00:00 2001 From: Chris Douglas <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 12:51:22 -0400 Subject: [PATCH 6/7] grammar Co-authored-by: Pierre Jolivet --- examples/douglas_jolivet_2026/README_structure.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/douglas_jolivet_2026/README_structure.md b/examples/douglas_jolivet_2026/README_structure.md index f28730d..788cab0 100644 --- a/examples/douglas_jolivet_2026/README_structure.md +++ b/examples/douglas_jolivet_2026/README_structure.md @@ -134,7 +134,7 @@ ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HyA2.mode - ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Hy_34.base -fo structure_HyA3 -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_HyA3.mode -fo structure_HyA3 -param P -tgv -2 -zero 1 ``` -- Bifurcations from $R_z$ to $A$ state: +- Bifurcations from $R_z$ to $A$ states: ```sh ff-mpirun -np $nproc modecompute.md -v 0 -dir $workdir -fi structure_Rz_23.base -fo structure_RzA -eps_target 0+0i -eps_nev 1 -strict 1 -eps_gen_hermitian -sym 1,1 ff-mpirun -np $nproc hopfcompute.md -v 0 -dir $workdir -fi structure_RzA.mode -fo structure_RzA -param P -tgv -2 -zero 1 From 277e8123bb21a7a3bd763fbb80166d14107af979 Mon Sep 17 00:00:00 2001 From: Chris Douglas <68232338+cmdoug@users.noreply.github.com> Date: Tue, 12 May 2026 12:51:30 -0400 Subject: [PATCH 7/7] typo Co-authored-by: Pierre Jolivet --- .../settings_douglas_jolivet_2026_structure.idp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp index fa9d799..f88658a 100644 --- a/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp +++ b/examples/douglas_jolivet_2026/settings_douglas_jolivet_2026_structure.idp @@ -37,7 +37,7 @@ // Define quantities for mesh adaptation and plotting in Paraview macro adaptu(u)[u, u#y, u#z]//EOM macro adaptf(f)[f, f#y, f#z]//EOM -// Name and order for real Paraview outputs +// Name and order for real ParaView outputs string ParaviewDataName = "disp"; string ParaviewDataNamef = ParaviewDataName; int[int] ParaviewOrder = [1];