Merge pull request #44 from UM-Bridge/cookies-problem

added cookies problem benchmark

benchmarks/cookies-problem/README.md

@@ -0,0 +1,73 @@
+# The cookies problem forward UQ benchmark
+
+## Overview
+
+This benchmark runs a forward uncertainty quantification problem for the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md) using the [Sparse Grids Matlab Kit]((https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit)) interface to UM-Bridge. See below for full description.
+
+## Authors
+- [Massimiliano Martinelli](mailto:[email protected])
+- [Lorenzo Tamellini](mailto:[email protected])
+
+## Run
+```
+docker run -it -p 4242:4242 linusseelinger/cookies-problem
+```
+
+## Properties
+
+Model     | Description
+---       | ---
+benchmark | sets the config options for the forward UQ benchmark (see below)
+
+### Benchmark configuration
+
+Mapping | Dimensions   | Description
+---     |---           |---
+input   | [8]          | These values modify the conductivity coefficient in the 8 cookies. They are i.i.d. uniform random variables in the range [-0.99 -0.2] (software does not check that inputs are within the bound) 
+output  | \[1\]        | The integral of the solution over the central subdomain (see definition of $$\Psi$$ at [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md) for info)
+
+Feature       | Supported
+---           |---
+Evaluate      | True
+Gradient      | False
+ApplyJacobian | False
+ApplyHessian  | False
+
+Config        | Type    | Default value   	| Can be changed in benchmark model	| Description
+---           |---      |---      		|---		  			| ---	
+NumThreads    | integer | 1     		| yes					| number of physical cores to be used by the solver 
+BasisDegree   | integer | 4       		| no					| Default degree of spline basis (must be a positive integer)
+Fidelity      | integer | 2       		| no					| Controls the number of mesh elements (must be a positive integer, see below for details)
+
+
+## Mount directories
+Mount directory | Purpose
+---             |---
+None            | 
+
+## Source code
+
+[Benchmark sources available at this folder.](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem)
+
+## Description
+
+![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/cookies-problem/cookies_domain.png "geometry of the cookies problem")
+
+The benchmark implements a forward uncertainty quantification problem for the [cookies model](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem/README.md). More specifically, we assume that the uncertain parameters $$y_n$$ appearing in the definition of the diffusion coefficient are uniform i.i.d. random variables on the range $$[-0.99, -0.2]$$ and we aim at computing the expected value of the quantity of interest (i.e., output of the model) $$\Psi$$, which is defined as the integral of the solution over $$F$$.
+
+The PDE is solved with an IGA solver that uses as basis splines of degree $$p=4$$ and maximal regularity, i.e. of continuity $$3$$, and the mesh has $$200 \times 200$$ elements (i.e., the fidelity config parameter is set to $$2$$). The structure of this benchmark is identical to the one discussed in \[1\]; however, raw numbers are different since in \[1\] the PDE solver employed was different (standard FEM with piecewise linear basis) and the mesh was also different.
+
+
+As a reference value, we provide the approximation of the expected value computed with a standard Smolyak sparse grid, based on Clenshaw--Curtis points, for level $$w=5$$, see e.g. \[2\]. The resulting sparse grid has 15713 points, and the corresponding approximation of the expected value is $$0.064196096847169$$.
+
+
+The script available [here](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem/run_forward_benchmark_in_matlab.m) generates the results, using the Sparse Grids Matlab Kit \[2\] for generating sparse grids. The Grids Matlab Kit is available on Github [here](https://github.com/lorenzo-tamellini/sparse-grids-matlab-kit) and a dedicated website with full resources including user manual is available [here](https://sites.google.com/view/sparse-grids-kit). 
+
+
+## Bibliography
+1 Joakim Bäck, Fabio Nobile, Lorenzo Tamellini, Raul Tempone, **Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: a numerical comparison**. In *Spectral and High Order Methods for Partial Differential Equations*, Vol. 76 of Lecture Notes in Computational Science and Engineering, Springer, 2011
+2 Chiara Piazzola, Lorenzo Tamellini, **The Sparse Grids Matlab Kit - a Matlab implementation of sparse grids for high-dimensional function approximation and uncertainty quantification**. ACM Transactions on Mathematical Software, 2023. 
+
+
+
+

benchmarks/cookies-problem/run_forward_benchmark_in_matlab.m

@@ -0,0 +1,123 @@
+%% ----------------------- forward UQ benchmark --------------------------
+% 
+% -------------------------------------------------------------------------
+
+
+%% preliminar operations:
+
+clear 
+
+% 1) run docker container with a command like
+% sudo docker run -it -p 4242:4242 <image name>
+
+
+% 2) add sparse grids matlab kit and um-bridge matlab client to path
+% addpath(genpath('~/GIT_projects/Github/sparse-grids-matlab-kit/')) % this is version 23.5 Robert
+% addpath(genpath('~/GIT_projects/Github/umbridge/matlab/'))
+
+
+% 3) to save results on file (sparse grids in .txt / .mat, results in .txt / .mat), set the flag below to true
+saving_stuff = true;
+
+
+% 4) to plot results, set the flag below to true
+plotting = true;
+
+
+%% setup model 
+
+% specify model port
+uri = 'http://0.0.0.0:4242';
+model = HTTPModel(uri,'benchmark');
+
+% config cookie solver with num threades
+config = struct('NumThreads',4);
+
+
+% wrap model in an @-function too
+Psi_fun = @(y) model.evaluate(y',config);
+
+% a simple call to test that things are working fine
+% y_test = [-0.5; -0.5; -0.5; -0.5; -0.5; -0.5; -0.5; -0.5;];
+% Psi_fun(y_test)
+
+%% create sparse grid
+
+% setup sparse grid ingredients. Instead of building the final grid in one go, we build it gradaully using a for
+% loop that creates a sequence of nested grids
+
+N = 8;
+knots = @(n) knots_CC(n,-0.99,-0.2);
+lev2knots = @lev2knots_doubling;
+idxset_rule = @(i) sum(i-1);
+idxset_level_max = 1; % <--- controls max size of sparse grid
+
+
+% To recycle evaluations from one grid to the next, we need some containers
+
+S_old         = []; % the previous sparse grid in extended format
+Sr_old        = []; % the previous sparse grid in reduced format
+Psi_evals_old = []; % the evaluations of the model output on the previous sparse grid
+
+% more containers, to save results at each iteration
+
+nb_pts = []; % nb of points of each sparse grid
+Psi_EV = []; % the expected value of Psi computed on each grid
+
+% here we go with the loop
+for w = 0:idxset_level_max
+
+    disp('========================================')
+    disp(strcat('w=',num2str(w)))
+    disp('========================================')
+
+    % create sparse grid in extended format. Recycle some work from previous sparse grid
+    S = create_sparse_grid(N,w,knots,lev2knots,idxset_rule,S_old);
+
+    % reduce sparse grid
+    Sr = reduce_sparse_grid(S);
+
+    % eval Psi on each point of sparse grid. Recycle available evaluations whenever possible
+    Psi_evals = evaluate_on_sparse_grid(Psi_fun,S,Sr,Psi_evals_old,S_old,Sr_old);
+
+    if saving_stuff
+        % save it to file. Type help export_sparse_grid_to_file for info on the saving format
+        grid_filename = strcat('sparse_grid_w=',num2str(w),'.txt');
+        export_sparse_grid_to_file(Sr,grid_filename,'with_weights');
+        % save also some minimal information in .mat format
+        grid_filename_mat = grid_filename(1:end-4); % i.e., remove .txt
+        save(grid_filename_mat,'S','Sr','Psi_evals')
+    end
+
+    % compute expected value and append the new value to previous container. Same for nb_pts
+    Psi_EV(end+1) = quadrature_on_sparse_grid(Psi_evals,Sr);    %#ok<SAGROW>
+    nb_pts(end+1) = Sr.size;                                    %#ok<SAGROW>
+
+    % update containers
+    S_old = S;
+    Sr_old = Sr;
+    Psi_evals_old = Psi_evals;
+
+end
+
+
+if saving_stuff
+    % save values to be plotted on file, in txt format. The data saved consists of the two vectors nb_pts and Psi_EV,
+    % they will be stored as rows of the txt file
+    save('cookies-benchmark-output.txt','nb_pts','Psi_EV','-ascii', '-double')
+    % save also on .mat
+    save('cookies-benchmark-output','nb_pts','Psi_EV')
+end
+
+
+if plotting
+    % in principle, only nb-pts and Psi_EV are needed, so if you have them on file already you can just load results 
+    % load('cookies-benchmark-output')
+    figure
+    semilogx(nb_pts,Psi_EV,'-ok','LineWidth',2,'MarkerFaceColor','k','DisplayName','sparse grid approx. of $\mathbf{E}[\Psi]$')
+    grid on
+    xlabel('sparse grids points')
+    ylabel('$\mathbf{E}[\Psi]$','interpreter','latex','rotation',0)
+    legend show
+    set(legend,'interpreter','latex','location','northwest','fontsize',14)
+end

models/cookies-problem/Dockerfile

@@ -0,0 +1,59 @@
+FROM ubuntu:mantic
+
+RUN apt update
+RUN DEBIAN_FRONTEND=noninteractive apt install -y bash python3-pip pbzip2 wget clang-17 git cmake curl zip pkg-config gfortran lld libomp5-17 libomp-17-dev
+RUN pip3 install --break-system-packages umbridge 
+
+
+RUN wget -O- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \ | gpg --dearmor | tee /usr/share/keyrings/oneapi-archive-keyring.gpg > /dev/null
+RUN echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" | tee /etc/apt/sources.list.d/oneAPI.list
+
+RUN apt update
+RUN apt install -y intel-basekit 
+# RUN . /opt/intel/oneapi/setvars.sh
+
+RUN git clone https://gitlab.com/max.martinelli/igatools.git && \
+  cd igatools && \
+  git submodule init && \
+  git submodule update && \
+  git checkout 576316cfcec2be3212cba7a95f2b2df60a3f64c2 && \
+  cd vcpkg && ./bootstrap-vcpkg.sh
+
+#  git checkout a335daa3bf480cb6e57e2ea8522e3817e841a66d && \
+
+
+WORKDIR /build_igatools
+RUN . /opt/intel/oneapi/setvars.sh && cd /build_igatools && \
+  cmake /igatools \
+  -G"Unix Makefiles" \
+  -DIGATOOLS_COMPONENT_DOCUMENTATION=OFF \
+  -DIGATOOLS_WITH_GMP=ON \
+  -DIGATOOLS_WITH_MKL=ON \
+  -DMKL_DIR=/opt/intel/oneapi/mkl/latest \
+  -DIGATOOLS_WITH_QUADRUPLE_PRECISION=OFF \
+  -DIGATOOLS_WITH_SERIALIZATION=OFF \
+  -DIGATOOLS_WITH_SUITESPARSE=OFF \
+  -DIGATOOLS_WITH_SUPERLU=OFF \
+  -DIGATOOLS_WITH_TRILINOS=OFF \
+  -DTBB_DIR=/opt/intel/oneapi/tbb/latest \
+  -DTBB_LIBRARY_DIR=/opt/intel/oneapi/tbb/latest/lib/intel64/gcc4.8/ \
+  -DTRILINOS_DIR=/opt/trilinos/current \
+  -DIGATOOLS_WITH_THREADS=ON \
+  -Wno-dev \
+  -DCMAKE_INSTALL_PREFIX=igatools \
+  -DCMAKE_BUILD_TYPE=Release \
+  -DIGATOOLS_STATIC_EXECUTABLE=OFF \
+  -DCMAKE_C_COMPILER=clang-17 \
+  -DCMAKE_CXX_COMPILER=clang++-17 \
+  -DBUILD_SHARED_LIBS=ON
+
+RUN make -j8 && \
+  make setup_tests && \
+  cd tests/ && \
+  make poisson_lorenzo.release 
+
+WORKDIR /
+
+COPY umbridge-server.py /
+
+CMD python3 umbridge-server.py

models/cookies-problem/README.md

@@ -0,0 +1,102 @@
+# The cookies model
+
+## Overview
+This model implements the so-called 'cookies problem' or 'cookies in the oven problem' \[1,2,3\], i.e., a simplified thermal equation in which the conductivity coefficient is uncertain in 8 circular subdomains ('the cookies'), whereas it is known (and constant) in the remaining of the domain ('the oven'). The PDE is solved by an isogeometric solver with maximum continuity splines, whose degree can be set by the user. See below for full description. 
+
+
+## Authors
+- [Massimiliano Martinelli](mailto:[email protected])
+- [Lorenzo Tamellini](mailto:[email protected])
+
+## Run
+```
+docker run -it -p 4242:4242 linusseelinger/cookies-problem
+```
+
+## Properties
+
+Model     | Description
+---       | ---
+forward   | forward evaluation of the cookies model, all config options can be modified by the user (see below)
+benchmark | sets the config options for the forward UQ benchmark [(see benchmark page)](https://github.com/UM-Bridge/benchmarks/tree/main/benchmarks/cookies-problem/README.md)
+
+### Forward
+
+Mapping | Dimensions | Description
+---     |---         |---
+input   | [8]        | These values modify the conductivity coefficient in the 8 cookies, each of them must be greater than -1 (software does not check that input values are valid)  
+output  | [1]        | The integral of the solution over the central subdomain (see definition of $$\Psi$$ below)
+
+Feature       | Supported
+---           |---
+Evaluate      | True
+Gradient      | False
+ApplyJacobian | False
+ApplyHessian  | False
+
+Config        | Type    | Default | Description
+---           |---      |---      |---
+NumThreads    | integer | 1       | number of physical cores to be used by the solver
+BasisDegree   | integer | 4       | Default degree of spline basis (must be a positive integer)
+Fidelity      | integer | 2       | Controls the number of mesh elements (must be a positive integer, see below for details)
+
+
+## Mount directories
+Mount directory | Purpose
+---             |---
+None            | 
+
+## Source code
+
+[Model sources here.](https://github.com/UM-Bridge/benchmarks/tree/main/models/cookies-problem)
+
+## Description
+
+![cookies-problem](https://raw.githubusercontent.com/UM-Bridge/benchmarks/main/models/cookies-problem/cookies_domain.png "geometry of the cookies problem")
+
+The model implements the version of the cookies problem in \[1\], see also e.g. \[2,3\] for slightly different versions. With reference to the computational domain $$D=[0,1]^2$$ in the figure above, the cookies model consists in the thermal diffusion problem below, where $$\mathbf{y}$$ are the uncertain parameters discussed in the following and $$\mathrm{x}$$ are physical coordinates 
+
+$$-\mathrm{div}\Big[ a(\mathbf{x},\mathbf{y}) \nabla u(\mathbf{x},\mathbf{y}) \Big] = f(\mathrm{x}), \quad \mathbf{x}\in D$$
+
+with homogeneous Dirichlet boundary conditions and forcing term defined as
+
+$$f(\mathrm{x}) = \begin{cases} 
+100 &\text{if } \,  \mathrm{x} \in F \\
+0 &\text{otherwise} 
+\end{cases}$$
+
+where $$F$$ is the square $$[0.4, 0.6]^2$$. The 8 subdomains with uncertain diffusion coefficient (the cookies) are circles with radius 0.13 and the following center coordinates:
+
+cookie | 1   | 2   | 3   | 4   | 5   | 6   | 7   | 8   |
+--     | --  | --  | --  | --  | --  | --  | --  | --  |
+x      | 0.2 | 0.5 | 0.8 | 0.2 | 0.8 | 0.2 | 0.5 | 0.8 |
+y      | 0.2 | 0.2 | 0.2 | 0.5 | 0.5 | 0.8 | 0.8 | 0.8 |
+
+The uncertain diffusion coefficient is defined as
+
+$$a = 1 + \sum_{i=1}^8 y_n \chi_n(\mathrm{x})$$ 
+
+where $$y_n>-1$$ and $$\chi_n(\mathrm{x}) = \begin{cases} 1 &\text{inside the n-th cookie} \\ 0 &\text{otherwise} \end{cases}$$
+
+
+The output of the simulation is the integral of the solution over $$F$$, i.e. $$\Psi = \int_F u(\mathrm{x}) d \mathrm{x}$$
+
+
+The PDE is solved with an IGA solver (see e.g. \[4\]) that uses as basis splines of degree $$p$$ (tunable by the user, default $$p=4$$) of maximal regularity, i.e. of continuity $$p-1$$. The computational mesh is an $$N\times N$$ quadrilateral mesh (cartesian product of knot lines) with square elements, with $$N=100 \times \mathrm{Fidelity}$$. The implementation is done using the C++ library IGATools \[5\], available at [gitlab.com/max.martinelli/igatools](gitlab.com/max.martinelli/igatools).  
+
+
+
+
+
+
+## Bibliography
+1 Joakim Bäck, Fabio Nobile, Lorenzo Tamellini, Raul Tempone, **Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: a numerical comparison**. In *Spectral and High Order Methods for Partial Differential Equations*, Vol. 76 of Lecture Notes in Computational Science and Engineering, Springer, 2011
+2 Jonas Ballani, Lars Grasedyck, **Hierarchical Tensor Approximation of Output Quantities of Parameter-Dependent PDEs**. *SIAM/ASA Journal of Uncertainty Quantification*, 2015
+3 Daniel Kressner, Christine Tobler, **Low-rank tensor Krylov subspace methods for parametrized linear systems**. *SIAM journal on matrix analysis and applications*, 2011
+4 Lourenco Beirao Da Veiga, Annalisa Buffa, Giancarlo Sangalli, Rafael Vazquez, **Mathematical analysis of variational isogeometric methods}**. *Acta Numerica*, 2014
+5 Miguel Sebastian Pauletti, Massimiliano Martinelli, Nicola Cavallini, Pablo Antolín, **IGATools: An isogeometric analysis library**. *SIAM journal on scientific computing*, 2015
+
+
+
+
+

models/cookies-problem/draw_geometry.m

@@ -0,0 +1,28 @@
+clear
+
+% basic domain
+
+% circles
+xc = [0.2 0.5 0.8 0.2 0.8 0.2 0.5 0.8];
+yc = [0.2 0.2 0.2 0.5 0.5 0.8 0.8 0.8];
+N = length(xc);
+
+radius = 0.13;
+theta = linspace(0,2*pi,200);
+
+for n=1:N
+    plot(xc(n) + radius*cos(theta), yc(n) + radius*sin(theta),'-k','LineWidth',2);
+    hold on
+    text(xc(n)-0.05,yc(n)+0.01,strcat('# ',num2str(n)),'FontSize',15)
+end
+
+% central square
+xsquare = [0.4 0.6 0.6 0.4 0.4];
+ysquare = [0.4 0.4 0.6 0.6 0.4];
+plot(xsquare,ysquare,'-k','LineWidth',2);
+text(0.485,0.5,'F','FontSize',15)
+
+axis square
+
+saveas(gcf,'cookies_domain','png')
+