# Benchmarking FBPINNs
---

This branch tests FBPINNs on various benchmarks. FBPINNs is proposed by the paper *[Finite Basis Physics-Informed Neural Networks (FBPINNs): a scalable domain decomposition approach for solving differential equations](https://arxiv.org/abs/2107.07871), B. Moseley, T. Nissen-Meyer and A. Markham, Jul 2021 ArXiv*. This code is extends the [repository implementing FBPINNs authored by bmoseley](https://github.com/benmoseley/FBPINNs) affiliated to the previous paper.

---

## Overview

This code extends the aforementioned codebase to support:

1. Irregular domains, by masking out unwanted regions from the original hyperrectangular domain
2. Both hard boundary constraint sand soft boundary constraint (boundary loss). When defining your own testcase, you must specify how to sample points from the boundary and how to calculate boundary loss.
3. Inverse problems, by supporting data loss.

The FBPINNs' loss is given by
$$
L_{FBPINNs} = L_{domain} + W_{B}L_{boundary}+W_{D}L_{data}
$$
where the domain loss is the physics loss (PDE residue squared), the boundary loss is the boundary conditions' equation residue squared, the data loss is the MSE between predicted and observed solutions at given points. If not defined by the specific PDE, the loss type will not be included in the total loss calculation.

## Environment

Please refer to the "installation" section of the aforementioned repository (FBPINNs by bmoselely).

## How to run this code on existing PDEs

### Configurations

Define the batch size, batch size (test) and subdomain division for each test case in `run_all_config.json`. The file already contains the default parameters.

You can configure other training hyperparameters in `constants.py`. The entry point of code `run_all.py` sets many parameters by creating and assigning values to a instance of the Constants class, but values not explicitly assigned in `run_all.py` take the default value as specified in `constants.py`. You might want to modify both files.

1. RANDOM (bool): whether to sample domain points randomly (if not, sample points on a uniform grid)
2. LRATE (float): learning rate
3. BOUNDARY_WEIGHT: weight for boundary loss, W_B
4. DATALOSS_WEIGHT: weight for data loss, W_D
5. SEED: random seed for torch and numpy
6. N_HIDDEN, N_LAYERS: net size for each subdomain's neural network 

### Run

To run a testcase whose name is  `[CaseName]`, run  `python runs/run_all.py [CaseName] ` in the current (fbpinns/) directory. The testcase is run three times with different seeds by default. For a list of casename see `run_all_config.json `

To run many testcases simultaneously (especially on a server with many CPUs), you can try to use the bash script `bash runs/run_all.sh`. In the future, we can modify the bash script, add "CUDA_VISIBLE_DEVICES=..." to run_all.sh to support running many testcases simultaneously on different GPUs.

### Output

After code finishes running:

1. Detailed output (generated by original code of bmoseley) can be found at `./results/`
2. A brief summary (plain text) can be found at `./benchmark_results/fb/`. The single-line file (filename contains the name of the testcase) corresponding to this run contains several space separated values. From first to last, they are:
   1. Time usage (seconds)
   2. Training loss
   3. L2 relative error
   4. L1 relative error
   5. Mean square error
   6. Mean absolute error
   7. Max error
   8. CSV error (as defined in the work "pdebench")
   9. fRMSE, which has three values: error in the low, mid, high (frequency)
   10. Subdomain division (of every dimension)
3. Visualization of the final solution (a png file) under `./benchmark_results/figs/fb/`
4. You can run `python runs/compare_3.py` to summarize results of all testcases, calculate the mean and standard deviation of metrics, and output a .csv file.

## How to run this code for the PINNacle paper

Main experiments and results in the PINNacle paper can be reproduced by running run  `python runs/run_all.py [CaseName] ` in the current (fbpinns/) directory. General hyperparameters are hard-coded into the file run_all.py, while hyperparameters that varies with testcases are listed in `runs/run_all_config.json`. You should use the given hyperparameter values in order to reproduce main results of FBPINNs.

To run the hyperparameter experiments, run `python runs/run_hyperparam.py [Casename] [Type] [Choice] ` in the current (fbpinns/) directory. The argument "Type" can be "width" meaning overlap ratio, or "div" meaning domain division; the argument "Choice" is the hyperparameter choice (usually 0, 1, 2, 3) meaning that we use the first, second, third or fourth chosen value for the hyperparameter.

To run the parametric PDE experiments, run `python runs/run_parametric.py [Casename] [Choice] ` in the current (fbpinns/) directory. The argument "Choice" can be 0,1,2 or 3 meaning that we use the first ... fourth chosen value as the parameter for the PDE.

## How to define your own PDEs

Please first read "defining your own PDEs" in the readme of FBPINNs by bmoselely. Like the original repo, each PDE is a class that inherits the class _Problem in `problems.py`. All files defining PDEs should be put in `pdes/`.

The following functions `get_gradients`, `physics_loss` , `boundary_condition` and `exact_solution` should be defined according to the instruction in the readme of FBPINNs by bmoselely. The following functions have additional requirements:

### \_\_init\_\_

Apart from self.d, the following variables must also be defined:

1. self.bbox (list of floats or ints): the bounding box of the domain
2. self.num_js (int): the number of partial derivative terms that are returned by `get_gradients` (or passed to `physics_loss`).

If the dimensionality of the exact solution is different from the dimensionality of network output y (e.g. exact solution only has several dimensions available, or problem is a inverse problem), then you should define self.exact_dim_select (type: slice, Python build-in type). For example, if this value is slice(1,2), then L2 relative error will be computed between the exact solution and y[1:2].

### boundary_condition

`boundary_condition` defines the **hard boundary constraints** to be applied to the NN ansatz, not the boundary loss. If your PDE does not use any hard boundary constraints, then do not overwrite this function and it will default to no hard boundary constraints.

### exact_solution

 `exact_solution` defines the exact / ground truth solution of the PDE.

If your PDE has a **analytical solution**, then I recommend you overwrite this function in your class and implement the analytical solution. The code will use it to compute test loss and final metrics including L2 relative error. 

If your PDE has only a **numerical reference solution**, then I recommend you not to overwrite this function. Put your reference solution (.dat file, format to be specified later) in under `../ref/` and add a line like  `self.load_ref_data("heat_complex", timepde=(0, 3))` in \_\_init\_\_.py. The code will automatically interpolate the reference solution to the test batch to compute test loss; the code will use the data points of the given reference solution to compute final metrics including L2 relative error.

### mask_x

If your PDE has **irregular domain**, then define a function with name "mask_x".

Input a tensor of shape (batch_size, nd), return a list or tensor (element type: Boolean) of shape (batch_size) representing whether each point (1, nd) in the input tensor belongs to the irregular domain.

### sample_bd, bd_loss

If you use **boundary loss to weakly enforce boundary conditions**, then define function "sample_bd" and "bd_loss".

Function "sample_bd" takes a int representing the number of boundary points to be sampled (N_bd) as input. It returns a list or numpy list or torch tensor of the sampled boundary points, on which boundary loss will be calculated.

Function "bd_loss" is similar to function "physics_loss" regarding their input and return values. The difference is that bd_loss calculates the boundary loss on the sampled boundary points instead of the physics loss.

### sample_data, data_loss

If you use data loss (have observed data) (e.g. inverse problems), then define "sample_data" and "data_loss".

Function "sample_data" takes no input and returns a list or numpy list or tensor of the input coordinates of the observed data points.

Function "data_loss" is similar to function "physics_loss" regarding their input and return values. "data_loss" should implement the MSE loss between the predicted and observed solutions.

