diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 2f0055ac..c927928f 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -19,7 +19,7 @@ jobs: .github/workflows/dependencies/documentation.sh echo "Installing python packages for docs..." python3 -m pip install --upgrade pip - python3 -m pip install sphinx sphinx_rtd_theme breathe sphinxcontrib.bibtex docutils + python3 -m pip install sphinx sphinx_rtd_theme breathe sphinxcontrib.bibtex docutils sphinx-copybutton sphinx-design - name: Install and Build run: | diff --git a/Docs/Readme.sphinx b/Docs/Readme.sphinx index 7ed0e5f4..30b15b59 100644 --- a/Docs/Readme.sphinx +++ b/Docs/Readme.sphinx @@ -23,6 +23,10 @@ If you have conda, you can install the necessary packages as follows: 2. Type "conda install sphinx_rtd_theme" + 3. Type "conda install sphinx-design" + + 4. Type "conda install -c conda-forge sphinx-copybutton" + If you don't have a conda Python installation, you get one by doing the following: 1. Go to https://conda.io/miniconda.html and download the Python 3.6 64-bit (bash installer), "Miniconda3-latest-Linux-x86_64.sh". Save this script to your hard drive. diff --git a/Docs/source/HeatEquation_UQ.rst b/Docs/source/HeatEquation_UQ.rst new file mode 100644 index 00000000..b69a3208 --- /dev/null +++ b/Docs/source/HeatEquation_UQ.rst @@ -0,0 +1,209 @@ +.. _tutorials_uq: + +.. _guided_pytuq_integration: + +.. _pytuq_quickstart: + +UQ with PyTUQ +=========== + +.. admonition:: **Time to Complete**: 30-60 minutes + :class: note + + **What you will learn**: + - Install PyTUQ + - Run an AMReX+PyTUQ to perform sensitivity analysis + +Overview +-------- + +AMReX simulations deliver high-fidelity, accurate results for complex physics problems, but the effect on simulation results due to uncertain inputs can require hundreds or thousands of simulations, making comprehensive analysis computationally prohibitive. + +This tutorial demonstrates how to improve efficiency without sacrificing accuracy by using polynomial chaos expansions to build fast surrogate models that identify which input parameters truly matter. + +PyTUQ (Python interface to the UQ Toolkit) provides specialized tools for surrogate construction and global sensitivity analysis, enabling rapid parameter space exploration and dimensionality reduction for scientific applications. + +We demonstrate how to integrate PyTUQ with your AMReX application through a practical workflow; the AMReX-based heat equation tutorial is equipped to perform sensitivity analysis. + +In these examples we model the heat equation + +.. math:: \frac{\partial\phi}{\partial t} = D\nabla^2 \phi, + +with initial condition + +.. math:: \phi_0 = 1 + A e^{-r^2 / (2V)}, + +where ``r`` is the distance from the center of the domain, and with uncertain parameters ``diffusion_coeff`` (:math:`D`), ``init_amplitude`` (:math:`A`), and ``init_variance`` (:math:`V`). +The outputs of interest are the maximum temperature, mean temperature, standard deviation of temperature, and the temperature at a specified point. + +.. toctree:: + :maxdepth: 1 + + HeatEquation_UQ_MathematicalDetails + +Located in ``amrex-tutorials/ExampleCodes/UQ/HeatEquation``, this example illustrates a complete forward UQ workflow from parameter sampling randomized input parameters to perform sensitivity analysis. +By understanding this example, you will have a basis for understanding how to adapt this workflow to your own AMReX application. + +More specifically, you can directly compare/diff ``amrex-tutorials/ExampleCodes/UQ/HeatEquation/main.cpp`` against the original heat equation tutorial ``amrex-tutorials/GuidedTutorials/HeatEquation_Simple/main.cpp`` to see exactly what source code changes are made to the AMReX application in this example. + +Installation +------------ + +We now describe the installation and workflow process on a local workstation. +First, install pytuq using this script (based on information provided in `pytuq/README.md `_): + +.. code-block:: bash + :caption: Pytuq installation script + + #!/bin/bash + + # 1. Clone repositories + git clone https://github.com/sandialabs/pytuq.git + + # 2. Create a conda environment with python (optional, you can add to an existing env) + conda create --name pytuq + conda activate pytuq + conda install python=3.11 + + # 3. Install PyTUQ and requirements + cd pytuq + pip install -r requirements.txt + pip install . + conda install dill + + # 4. Verify installation + conda list | grep pytuq # Should show pytuq 1.0.0 + +Examples +-------- + +C++ AMReX + PyTUQ (BASH driven) +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + **Prerequisites**: + + - AMReX and pytuq cloned at the same directory level as amrex-tutorials + - pytuq installation described above + - GNU Parallel for task management: ``sudo apt-get install parallel`` (mpiexec option also exists in script) + + .. code-block:: bash + :caption: Build AMReX appplication + + cd /path/to/amrex-tutorials/ExampleCodes/UQ/HeatEquation/ + make -j4 + + .. code-block:: bash + :caption: Run sensitivity analysis with bash script + + ./workflow_uqpc.x + +Understanding GNU Parallel Workflow Pattern +------------------------------------------- + + The ``workflow_uqpc.x`` bash script relies on the user augmenting their codes to write outputs of interest to ASCII text files. + In this case, the ``main.cpp`` was modified from the ``amrex-tutorials/GuidedTutorials/HeatEquation_Simple/main.cpp`` in the following ways: + + First, support for parsing ``diffusion_coeff``, ``init_amplitude``, and ``init_variance`` from the input file and command line were added. + + Second, support for writing outputs of interest to ASCII text files is added. The ``datalog_filename`` input parameter is generated by the bash + script to give each simulation output a unique identifier. + The ``datalog_int`` parameter gives the user the option to write the outputs of interest at a given time step interval, but in this example, + the outputs of interest after the final step are those that matter, and are extracted by the bash script to create a master output file + containing a separate set of simulation outputs of interest in each row. + + The bash script calls PyTUQ scripts that generate an input parameter files for training and testing (``ptrain.txt`` and ``ptest.txt``) + based on polynomial chaos settings, and then uses GNU Parallel to run multiple simulations + efficiently and collect outputs into results files (``ytrain.txt`` and ``ytest.txt``) that PyTUQ can use for surrogate model fitting. + +Understanding the Output +------------------------ + +| `pnames.txt`, `outnames.txt` +| -names of input and output parameters + +| `param_margpc.txt` +| -Each row contains the mean and standard deviation for an uncertain input parameter + +| `qtrain.txt`, `qtest.txt` +| -each row is a separate set of normal random variables used to generate uncertain inputs + +| `ptrain.txt`, `ptest.txt`, `pall.txt` +| -each row is a separate set of input parameters for each simulation + +| `stdoutlog_train##.txt`, `stdoutlog_test##.txt` +| -All standard output from AMReX simulations for training and testing data. The numbers refer to separate simulations. + +| `datalog_train##.txt`, `datalog_test##.txt` +| -Specifically-chosen output from AMReX simulations for training and testing data. In this example it reports the outputs (max, mean, standard deviation, and specific cell temperature) at user-specified intervals. + +| `ytrain.txt`, `ytest.txt`, `yall.txt` +| -agglomeration of outputs of interest from all simulations + +| `results.pk` +| -Python pickle file encapsulating the results + +| `labels.txt` +| -list of labels of simulation types (training or testing) for diagnostic/plot generation purposes + +| `xx__.png` +| -scatter plot of 2 inputs for training and testing + +| `pcoord_1.png`, `pcoord_1_lab_Testing.png`, `pcoord_1_lab_Training.png` +| -graphical representation of how inputs correlate to outputs for each individual simulation + +| `yx_.png`, `yx__log.png` +| -scatter plots of output as a function of each input + +| `yxx_#.png` +| -scatter plots of output a function of multiple inputs + +| `pdf_tri_inputs.png`, `pdf_tri_output.png` +| -PDFs of inputs and outputs + +| `ensemble.png` +| -graphical display of all output values + +| `idm_#_training.png`, `idm_#_testing.png` +| -graphical display of output values + +| `dm_#.png` +| -”data vs model” parity plots for each output; compares the predicted values from a surrogate model or approximation against the true or actual values from the full computational model + +| `fit_s#_training.png`, `fit_s#_testing.png` +| -Shows model vs PC approximation for a single simulation. + +| `pdf_output_#.png`, `pdf_joyplot.png` +| -PDF of output variables + +| `allsens_main.txt`, `sens_main.png` +| -raw data and plot for parameter sensitivities + +| `jsens_#.png`, `Jsens_ave.png` +| -joint sensitivities of output with respect to inputs + +| `sensmat_main.png` +| -sensitivity matrix + +| `pcslices_o#.png` +| -polynomial chaos fits + +| `pccont_o#_d#_d#.png` +| -polynomial chaos fits of output variables with respect to two input variables + + +Additional Resources +-------------------- + +**PyTUQ Resources:** + +- `PyTUQ Documentation `_ +- `PyTUQ Examples directory `_ +- eebaill, ksargsyan, & Bert Debusschere. (2025). sandialabs/pytuq: v1.0.0z (v1.0.0z). Zenodo. https://doi.org/10.5281/zenodo.17110054 + +**AMReX Resources:** + +- `AMReX Documentation `_ + +**Uncertainty Quantification Theory:** + +- Ghanem, Roger, David Higdon, and Houman Owhadi, eds. *Handbook of Uncertainty Quantification*. Vol. 6. New York: Springer, 2017. (For workflow, plotting, and analysis specifics) diff --git a/Docs/source/HeatEquation_UQ_MathematicalDetails.rst b/Docs/source/HeatEquation_UQ_MathematicalDetails.rst new file mode 100644 index 00000000..6db430cb --- /dev/null +++ b/Docs/source/HeatEquation_UQ_MathematicalDetails.rst @@ -0,0 +1,48 @@ +.. _pytuq_mathematical_details: + +Parameter Sensitivity via Polynomial Chaos +------------------------------------------- + +In this example, PyTUQ constructs polynomial chaos expansions to analyze how uncertain physical parameters affect simulation outputs in a heat diffusion problem. + +The Heat Equation +^^^^^^^^^^^^^^^^^ + +The governing equation for this tutorial is the heat diffusion equation: + +.. math:: \frac{\partial\phi}{\partial t} = D\nabla^2 \phi, + +with initial condition + +.. math:: \phi_0 = 1 + A e^{-r^2 / (2V)}, + +where ``r`` is the distance from the center of the domain, and with uncertain parameters ``diffusion_coeff`` (:math:`D`), ``init_amplitude`` (:math:`A`), and ``init_variance`` (:math:`V`). + + + +Quantities of Interest (Outputs) +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The simulation extracts four statistical quantities at the final timestep: + +1. **max_temp**: Maximum temperature in the domain +2. **mean_temp**: Average temperature across all cells +3. **std_temp**: Standard deviation of temperature +4. **cell_temp**: Temperature in a particular cell, in this case, cell (9,9,9) + +These outputs are computed in ``main.cpp``and written to the datalog file. + +PyTUQ Workflow +^^^^^^^^^^^^^^ + +PyTUQ uses polynomial chaos expansion to construct a surrogate model: + +1. **Parameter sampling**: Generate sample points in the 3D parameter space based on the specified distributions +2. **Forward simulations**: Run the AMReX heat equation solver for each parameter set +3. **Surrogate construction**: Fit polynomial chaos coefficients to map inputs → outputs +4. **Sensitivity analysis**: Compute Sobol indices to identify which parameters most influence each output + +The connection is: + +- **Inputs**: ParmParse parameters (``diffusion_coeff``, ``init_amplitude``, ``init_variance``) specified in ``inputs`` file or command line +- **Outputs**: Quantities of interest extracted from datalog files diff --git a/Docs/source/conf.py b/Docs/source/conf.py index a8f544f5..f8519ad7 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -39,6 +39,8 @@ def get_amrex_version(): extensions = ['sphinx.ext.mathjax', 'sphinx.ext.githubpages', 'sphinx.ext.viewcode', + 'sphinx_design', + 'sphinx_copybutton', 'sphinx.ext.intersphinx', 'sphinx_rtd_theme'] @@ -51,6 +53,9 @@ def get_amrex_version(): # Add any paths that contain templates here, relative to this directory. templates_path = ['ytemplates'] +# sphinx-copybutton configuration +copybutton_exclude = '.linenos, .gp, .go' + # The suffix(es) of source filenames. # You can specify multiple suffix as a list of string: # diff --git a/Docs/source/index.rst b/Docs/source/index.rst index f185fc6a..cba245ae 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -50,6 +50,7 @@ sorted by the following categories: - :ref:`GPU` -- Offload work to the GPUs using AMReX tools. - :ref:`Linear Solvers` -- Examples of several linear solvers. - :ref:`ML/PYTORCH` -- Use of pytorch models to replace point-wise computational kernels. +- :ref:`UQ with pytuq` -- Use of pytuq methods to perform generalized sensitivity analysis. - :ref:`MPMD` -- Usage of AMReX-MPMD (Multiple Program Multiple Data) framework. - :ref:`MUI` -- Incorporates the MxUI/MUI (Multiscale Universal interface) frame into AMReX. - :ref:`Particles` -- Basic usage of AMReX's particle data structures. @@ -72,6 +73,7 @@ sorted by the following categories: GPU_Tutorial LinearSolvers_Tutorial ML_Tutorial + HeatEquation_UQ MPMD_Tutorials MUI_Tutorial Particles_Tutorial @@ -97,6 +99,8 @@ sorted by the following categories: .. _`Linear Solvers`: LinearSolvers_Tutorial.html +.. _`UQ`: UQ_Tutorial.html + .. _`MPMD`: MPMD_Tutorials.html .. _`MUI`: MUI_Tutorial.html diff --git a/ExampleCodes/UQ/HeatEquation/GNUmakefile b/ExampleCodes/UQ/HeatEquation/GNUmakefile new file mode 100644 index 00000000..90bcb595 --- /dev/null +++ b/ExampleCodes/UQ/HeatEquation/GNUmakefile @@ -0,0 +1,16 @@ +# AMREX_HOME defines the directory in which we will find all the AMReX code. +AMREX_HOME ?= ../../../../amrex + +DEBUG = FALSE +USE_MPI = TRUE +USE_OMP = FALSE +COMP = gnu +DIM = 3 + +include $(AMREX_HOME)/Tools/GNUMake/Make.defs + +include ./Make.package + +include $(AMREX_HOME)/Src/Base/Make.package + +include $(AMREX_HOME)/Tools/GNUMake/Make.rules diff --git a/ExampleCodes/UQ/HeatEquation/Make.package b/ExampleCodes/UQ/HeatEquation/Make.package new file mode 100644 index 00000000..7f43e5e8 --- /dev/null +++ b/ExampleCodes/UQ/HeatEquation/Make.package @@ -0,0 +1,2 @@ +CEXE_sources += main.cpp + diff --git a/ExampleCodes/UQ/HeatEquation/inputs b/ExampleCodes/UQ/HeatEquation/inputs new file mode 100644 index 00000000..2a654017 --- /dev/null +++ b/ExampleCodes/UQ/HeatEquation/inputs @@ -0,0 +1,16 @@ +# Grid/Domain parameters +n_cell = 32 # number of cells on each side of domain +max_grid_size = 16 # size of each box (or grid) + +# Time stepping parameters +nsteps = 100 # total steps in simulation +dt = 5.e-8 # time step + +# Output control +plot_int = -1 # how often to write plotfile (-1 = no plots) +datalog_int = -1 # how often to write datalog (-1 = no regular output) + +# Physics parameters (these are what we vary for UQ and values here will be overwritten +diffusion_coeff = 100. # diffusion coefficient for heat equation +init_amplitude = 1. # amplitude of initial temperature profile +init_variance = 0.01 # variance of initial temperature profile diff --git a/ExampleCodes/UQ/HeatEquation/main.cpp b/ExampleCodes/UQ/HeatEquation/main.cpp new file mode 100644 index 00000000..1f383936 --- /dev/null +++ b/ExampleCodes/UQ/HeatEquation/main.cpp @@ -0,0 +1,351 @@ +/* + * A simplified single file version of the HeatEquation_EX0_C exmaple. + * This code is designed to be used with Demo_Tutorial.rst. + * + */ + + +#include +#include +#include + + +int main (int argc, char* argv[]) +{ + amrex::Initialize(argc,argv); + { + + // ********************************** + // DECLARE SIMULATION PARAMETERS + // ********************************** + + // number of cells on each side of the domain + int n_cell; + + // size of each box (or grid) + int max_grid_size; + + // total steps in simulation + int nsteps; + + // how often to write a plotfile + int plot_int; + + // time step + amrex::Real dt; + + // ********************************** + // DECLARE PHYSICS PARAMETERS + // ********************************** + + // diffusion coefficient for heat equation + amrex::Real diffusion_coeff; + + // amplitude of initial temperature profile + amrex::Real init_amplitude; + + // variance of the intitial temperature profile + amrex::Real init_variance; + + // ********************************** + // DECLARE DATALOG PARAMETERS + // ********************************** + const int datwidth = 24; + const int datprecision = 16; + const int timeprecision = 13; + int datalog_int = -1; // Interval for regular output (<=0 means no regular output) + std::string datalog_filename = "datalog.txt"; + + // ********************************** + // READ PARAMETER VALUES FROM INPUT DATA + // ********************************** + // inputs parameters + { + // ParmParse is way of reading inputs from the inputs file + // pp.get means we require the inputs file to have it + // pp.query means we optionally need the inputs file to have it - but we must supply a default here + amrex::ParmParse pp; + + // We need to get n_cell from the inputs file - this is the number of cells on each side of + // a square (or cubic) domain. + pp.get("n_cell",n_cell); + + // The domain is broken into boxes of size max_grid_size + pp.get("max_grid_size",max_grid_size); + + // Default nsteps to 10, allow us to set it to something else in the inputs file + nsteps = 10; + pp.query("nsteps",nsteps); + + // Default plot_int to -1, allow us to set it to something else in the inputs file + // If plot_int < 0 then no plot files will be written + plot_int = -1; + pp.query("plot_int",plot_int); + + // time step + pp.get("dt",dt); + + // Default datalog_int to -1, allow us to set it to something else in the inputs file + // If datalog_int < 0 then no plot files will be written + datalog_int = -1; + pp.query("datalog_int",datalog_int); + + datalog_filename = "datalog.txt"; + pp.query("datalog",datalog_filename); + + // ********************************** + // READ PHYSICS PARAMETERS + // ********************************** + + // Diffusion coefficient - controls how fast heat spreads + diffusion_coeff = 1.0; + pp.query("diffusion_coeff", diffusion_coeff); // Note: input name is "diffusion" but variable is "diffusion_coeff" + + // Initial temperature amplitude + init_amplitude = 1.0; + pp.query("init_amplitude", init_amplitude); + + // Initial temperature variance + // Smaller values = more concentrated, larger values = more spread out + init_variance = 0.01; + pp.query("init_variance", init_variance); + } + + // ********************************** + // DEFINE SIMULATION SETUP AND GEOMETRY + // ********************************** + + // make BoxArray and Geometry + // ba will contain a list of boxes that cover the domain + // geom contains information such as the physical domain size, + // number of points in the domain, and periodicity + amrex::BoxArray ba; + amrex::Geometry geom; + + // define lower and upper indices + amrex::IntVect dom_lo(0,0,0); + amrex::IntVect dom_hi(n_cell-1, n_cell-1, n_cell-1); + + // Make a single box that is the entire domain + amrex::Box domain(dom_lo, dom_hi); + + // Initialize the boxarray "ba" from the single box "domain" + ba.define(domain); + + // Break up boxarray "ba" into chunks no larger than "max_grid_size" along a direction + ba.maxSize(max_grid_size); + + // This defines the physical box, [0,1] in each direction. + amrex::RealBox real_box({ 0., 0., 0.}, + { 1., 1., 1.}); + + // periodic in all direction + amrex::Array is_periodic{1,1,1}; + + // This defines a Geometry object + geom.define(domain, real_box, amrex::CoordSys::cartesian, is_periodic); + + // extract dx from the geometry object + amrex::GpuArray dx = geom.CellSizeArray(); + + // Nghost = number of ghost cells for each array + int Nghost = 1; + + // Ncomp = number of components for each array + int Ncomp = 1; + + // How Boxes are distrubuted among MPI processes + amrex::DistributionMapping dm(ba); + + // we allocate two phi multifabs; one will store the old state, the other the new. + amrex::MultiFab phi_old(ba, dm, Ncomp, Nghost); + amrex::MultiFab phi_new(ba, dm, Ncomp, Nghost); + amrex::MultiFab phi_tmp(ba, dm, Ncomp, Nghost); // temporary used to calculate standard deviation + + // time = starting time in the simulation + amrex::Real time = 0.0; + + // ********************************** + // INITIALIZE DATA LOOP + // ********************************** + + // loop over boxes + for (amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi) + { + const amrex::Box& bx = mfi.validbox(); + + const amrex::Array4& phiOld = phi_old.array(mfi); + + // ********************************** + // SET INITIAL TEMPERATURE PROFILE + // ********************************** + // Formula: T = 1 + amplitude * exp(-r^2 / (2*variance)) + // where r is distance from center (0.5, 0.5, 0.5) + // + // Parameters: + // - amplitude: controls peak temperature above baseline (1.0) + // - variance: controls spread of initial hot spot + // - smaller variance = more concentrated + // - larger variance = more spread out + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) + { + + // ********************************** + // SET VALUES FOR EACH CELL + // ********************************** + + amrex::Real x = (i+0.5) * dx[0]; + amrex::Real y = (j+0.5) * dx[1]; + amrex::Real z = (k+0.5) * dx[2]; + amrex::Real rsquared = ((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)+(z-0.5)*(z-0.5))/(2*init_variance); + phiOld(i,j,k) = 1. + init_amplitude * std::exp(-rsquared); + }); + } + + // ********************************** + // WRITE DATALOG FILE + // ********************************** + if (amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream datalog(datalog_filename); // truncate mode to start fresh + datalog << "#" << std::setw(datwidth-1) << " max_temp"; + datalog << std::setw(datwidth) << " mean_temp"; + datalog << std::setw(datwidth) << " std_temp"; + datalog << std::setw(datwidth) << " cell_temp"; + datalog << std::endl; + datalog.close(); + } + + // ********************************** + // WRITE INITIAL PLOT FILE + // ********************************** + + // Write a plotfile of the initial data if plot_int > 0 + if (plot_int > 0) + { + int step = 0; + const std::string& pltfile = amrex::Concatenate("plt",step,5); + WriteSingleLevelPlotfile(pltfile, phi_old, {"phi"}, geom, time, 0); + } + + + // ********************************** + // MAIN TIME EVOLUTION LOOP + // ********************************** + + for (int step = 1; step <= nsteps; ++step) + { + // fill periodic ghost cells + phi_old.FillBoundary(geom.periodicity()); + + // new_phi = old_phi + dt * Laplacian(old_phi) + // loop over boxes + for ( amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi ) + { + const amrex::Box& bx = mfi.validbox(); + + const amrex::Array4& phiOld = phi_old.array(mfi); + const amrex::Array4& phiNew = phi_new.array(mfi); + + // advance the data by dt + amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE (int i, int j, int k) + { + + // ********************************** + // EVOLVE VALUES FOR EACH CELL + // ********************************** + + phiNew(i,j,k) = phiOld(i,j,k) + dt * diffusion_coeff * + ( (phiOld(i+1,j,k) - 2.*phiOld(i,j,k) + phiOld(i-1,j,k)) / (dx[0]*dx[0]) + +(phiOld(i,j+1,k) - 2.*phiOld(i,j,k) + phiOld(i,j-1,k)) / (dx[1]*dx[1]) + +(phiOld(i,j,k+1) - 2.*phiOld(i,j,k) + phiOld(i,j,k-1)) / (dx[2]*dx[2]) + ); + }); + } + + // find the value in cell (9,9,9) + amrex::ReduceOps reduce_op; + amrex::ReduceData reduce_data(reduce_op); + using ReduceTuple = typename decltype(reduce_data)::Type; + + for ( amrex::MFIter mfi(phi_old); mfi.isValid(); ++mfi ) + { + const amrex::Box& bx = mfi.validbox(); + + const amrex::Array4& phiOld = phi_old.array(mfi); + const amrex::Array4& phiNew = phi_new.array(mfi); + + // advance the data by dt + reduce_op.eval(bx, reduce_data, + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple + { + if (i==9 && j==9 && k==9) { + return{phiNew(i,j,k)}; + } else { + return {0.}; + } + }); + } + + amrex::Real cell_temperature = amrex::get<0>(reduce_data.value()); + amrex::ParallelDescriptor::ReduceRealSum(cell_temperature); + + // ********************************** + // INCREMENT + // ********************************** + + // update time + time = time + dt; + + // copy new solution into old solution + amrex::MultiFab::Copy(phi_old, phi_new, 0, 0, 1, 0); + + // Tell the I/O Processor to write out which step we're doing + amrex::Print() << "Advanced step " << step << "\n"; + + // ********************************** + // WRITE DATALOG AT GIVEN INTERVAL + // ********************************** + + // Check if we should write datalog + bool write_datalog = false; + if (step == nsteps) { + write_datalog = true; // Write final step + } else if (datalog_int > 0 && step % datalog_int == 0) { + write_datalog = true; // Write every datalog_int steps + } + + amrex::MultiFab::Copy(phi_tmp, phi_new, 0, 0, 1, 0); + amrex::Real max_temperature = phi_new.max(0); + amrex::Real mean_temperature = phi_new.sum(0) / phi_new.boxArray().numPts(); + phi_tmp.plus(-mean_temperature,0,1,0); + amrex::Real std_temperature = phi_tmp.norm2(0); // compute sqrt( sum(phi_tmp_i^2) ); + + if (write_datalog && amrex::ParallelDescriptor::IOProcessor()) { + std::ofstream datalog(datalog_filename, std::ios::app); + + // Write 4 statistics + datalog << std::setw(datwidth) << std::setprecision(datprecision) << max_temperature; + datalog << std::setw(datwidth) << std::setprecision(datprecision) << mean_temperature; + datalog << std::setw(datwidth) << std::setprecision(datprecision) << std_temperature; + datalog << std::setw(datwidth) << std::setprecision(datprecision) << cell_temperature; + datalog << std::endl; + + datalog.close(); + } + + // ********************************** + // WRITE PLOTFILE AT GIVEN INTERVAL + // ********************************** + + // Write a plotfile of the current data (plot_int was defined in the inputs file) + if (plot_int > 0 && step%plot_int == 0) + { + const std::string& pltfile = amrex::Concatenate("plt",step,5); + WriteSingleLevelPlotfile(pltfile, phi_new, {"phi"}, geom, time, step); + } + } + + } + amrex::Finalize(); + return 0; +} diff --git a/ExampleCodes/UQ/HeatEquation/workflow_uqpc.x b/ExampleCodes/UQ/HeatEquation/workflow_uqpc.x new file mode 100755 index 00000000..d70661fb --- /dev/null +++ b/ExampleCodes/UQ/HeatEquation/workflow_uqpc.x @@ -0,0 +1,166 @@ +#!/bin/bash -e +#===================================================================================== + +# indicate location of pytuq repo here +export PYTUQ_HOME=../../../../pytuq + +# relative location of pytuq and uqpc apps +export PUQAPPS=$PYTUQ_HOME/apps +export UQPC=$PYTUQ_HOME/apps/uqpc + +################################ +## 0. Setup the problem ## +################################ + +## Given mean and standard deviation of each normal random parameter +# inputs are diffusion_coeff, init_amplitude, init_variance +echo "400 100 " > param_margpc.txt +echo "1.0 0.25" >> param_margpc.txt +echo "0.01 0.0025" >> param_margpc.txt +# Hermite-Gaussian PC +PC_TYPE=HG +INPC_ORDER=1 +# Creates input PC coefficient file pcf.txt (will have lots of zeros since we assume independent inputs) +${PUQAPPS}/pc_prep.py -f marg -i param_margpc.txt -p ${INPC_ORDER} + +# Number of samples requested +NTRN=99 # Training +NTST=20 # Testing + +# Extract dimensionality d (i.e. number of input parameters) +DIM=`awk 'NR==1{print NF}' pcf.txt` + +# Output PC order +OUTPC_ORDER=3 + +############################### +## 1. Prepare the inputs ## +############################### + +# Prepare inputs for the black-box model (use input PC to generate input samples for the model) +# This creates files ptrain.txt, ptest.txt (train/test parameter inputs), qtrain.txt, qtest.txt (corresponding train/test stochastic PC inputs) +${PUQAPPS}/pc_sam.py -f pcf.txt -t ${PC_TYPE} -n $NTRN +mv psam.txt ptrain.txt; mv qsam.txt qtrain.txt +${PUQAPPS}/pc_sam.py -f pcf.txt -t ${PC_TYPE} -n $NTST +mv psam.txt ptest.txt; mv qsam.txt qtest.txt + +# create pnames.txt and outnames.txt with input parameter names and output QoI names +echo "diffusion_coef" > pnames.txt +echo "init_amplitude" >> pnames.txt +echo "init_variance" >> pnames.txt +echo "max_temp" > outnames.txt +echo "mean_temp" >> outnames.txt +echo "std_temp" >> outnames.txt +echo "cell_temp" >> outnames.txt + +################################ +## 2. Run the black-box model ## +################################ + +# Run the black-box model, both training and testing +# use parallel to launch multiple jobs at a time, each using 1 MPI rank +# ptrain.txt is N x d input matrix, each row is a parameter vector of size d +# ytrain.txt is N x o output matrix, each row is a output vector of size o +# Similar for testing +parallel --jobs 4 --keep-order --colsep ' ' './main3d.gnu.MPI.ex inputs diffusion_coeff={1} init_amplitude={2} init_variance={3} datalog=datalog_train{#}.txt > stdoutlog_train{#}.txt 2>&1 ; tail -1 datalog_train{#}.txt' :::: ptrain.txt > ytrain.txt +parallel --jobs 4 --keep-order --colsep ' ' './main3d.gnu.MPI.ex inputs diffusion_coeff={1} init_amplitude={2} init_variance={3} datalog=datalog_test{#}.txt > stdoutlog_test{#}.txt 2>&1 ; tail -1 datalog_test{#}.txt' :::: ptest.txt > ytest.txt + +################################ +# 2. Run the black-box model (alternative using mpiexec) +################################ + +## Initialize job index +#job_index=0 +# +## Read each line from ptrain.txt +#while IFS=' ' read -r diffusion_coeff init_amplitude init_variance; do +# +# # Prepare the datalog and stdout log filenames +# datalog_filename="datalog_train${job_index}.txt" +# stdoutlog_filename="stdoutlog_train${job_index}.txt" +# +# # Run the job with mpiexec +# mpiexec -n 4 ./main3d.gnu.MPI.ex inputs diffusion_coeff=$diffusion_coeff init_amplitude=$init_amplitude init_variance=$init_variance datalog=$datalog_filename $stdoutlog_filename 2>&1 +# +# # If successful, output the last line of the datalog +# if [ $? -eq 0 ]; then +# tail -1 "$datalog_filename" >> ytrain.txt +# fi +# job_index=$((job_index+1)) +#done < ptrain.txt +# +## Initialize job index +#job_index=0 +# +## Read each line from ptest.txt +#while IFS=' ' read -r diffusion_coeff init_amplitude init_variance; do +# +# # Prepare the datalog and stdout log filenames +# datalog_filename="datalog_test${job_index}.txt" +# stdoutlog_filename="stdoutlog_test${job_index}.txt" +# +# # Run the job with mpiexec +# mpiexec -n 4 ./main3d.gnu.MPI.ex inputs diffusion_coeff=$diffusion_coeff init_amplitude=$init_amplitude init_variance=$init_variance datalog=$datalog_filename $stdoutlog_filename 2>&1 +# +# # If successful, output the last line of the datalog +# if [ $? -eq 0 ]; then +# tail -1 "$datalog_filename" >> ytest.txt +# fi +# job_index=$((job_index+1)) +#done < ptest.txt + +############################## +# 3. Build PC surrogate ## +############################## + +# Build surrogate for each output (in other words, build output PC) +# This creates files results.pk (Python pickle file encapsulating the results) +${UQPC}/uq_pc.py -r offline -c pcf.txt -x ${PC_TYPE} -d $DIM -o ${INPC_ORDER} -m anl -s rand -n $NTRN -v $NTST -t ${OUTPC_ORDER} + +################################### +## 4. Visualize the i/o data ## +################################### + +awk '{print "Training"}' ytrain.txt > labels.txt +awk '{print "Testing"}' ytest.txt >> labels.txt +cat ytrain.txt ytest.txt > yall.txt +cat ptrain.txt ptest.txt > pall.txt + +${PUQAPPS}/plot_xx.py -x pall.txt -l labels.txt +${PUQAPPS}/plot_pcoord.py -x pall.txt -y yall.txt -l labels.txt + +${PUQAPPS}/plot_yx.py -x ptrain.txt -y ytrain.txt -c 3 -r 1 +${PUQAPPS}/plot_yxx.py -x ptrain.txt -y ytrain.txt +${PUQAPPS}/plot_pdfs.py -p ptrain.txt; cp pdf_tri.png pdf_tri_inputs.png +${PUQAPPS}/plot_pdfs.py -p ytrain.txt; cp pdf_tri.png pdf_tri_outputs.png +${PUQAPPS}/plot_ens.py -y ytrain.txt + +# A lot of .png files are created for visualizing the input/output data + +################################ +## 5. Postprocess the results ## +################################ + +# Plot model vs PC for each output +${UQPC}/plot.py dm training testing +# Plot model vs PC for each sample +${UQPC}/plot.py fit training testing + +# Plot output pdfs +${UQPC}/plot.py pdf +# Plot output pdfs in joyplot format +${UQPC}/plot.py joy + +# Plot main Sobol sensitivities (barplot) +${UQPC}/plot.py sens main +# Plot joint Sobol sensitivities (circular plot) +${UQPC}/plot.py jsens +# Plot matrix of Sobol sensitivities (rectangular plot) +${UQPC}/plot.py sensmat main + +# Plot 1d slices of the PC surrogate +${UQPC}/plot.py 1d training testing +# Plot 2d slices of the PC surrogate +${UQPC}/plot.py 2d + +# This creates .png files for visualizing PC surrogate results and sensitivity analysis diff --git a/GuidedTutorials/HeatEquation_Simple/main.cpp b/GuidedTutorials/HeatEquation_Simple/main.cpp index b3772626..b1f91aac 100644 --- a/GuidedTutorials/HeatEquation_Simple/main.cpp +++ b/GuidedTutorials/HeatEquation_Simple/main.cpp @@ -90,7 +90,7 @@ int main (int argc, char* argv[]) // This defines the physical box, [0,1] in each direction. amrex::RealBox real_box({ 0., 0., 0.}, - { 1., 1., 1.}); + { 1., 1., 1.}); // periodic in all direction amrex::Array is_periodic{1,1,1}; @@ -217,10 +217,7 @@ int main (int argc, char* argv[]) } } - } amrex::Finalize(); return 0; } - -