# set

The `set` command is used to set some properties of the analysis.

## Linear Elastic Simplification

In a general non-linear context, each substep requires at least two system solving, one for initial guess, one for\
convergence test.\
There is no general way to detect a linear system automatically,\
with which we know the initial guess leads to convergence so that further system solving is not necessary.\
To speed up simulation, one can use the following command to indicate the system is linear.

```
set linear_system true
```

With this enabled, only one system solving will be performed in each substep.

## Substepping Control

To use a fixed substep size, users can define

```
set fixed_step_size true
```

Otherwise, the algorithm will automatically substep the current step if convergence is not met. The time step control\
strategy cannot be customized.

To define the initial step size, users can use

```
set ini_step_size (1)
# (1) double, substep size
```

To define the maximum/minimum step size, users can use

```
set max_step_size (1)
set min_step_size (1)
# (1) double, target time step size
```

## Solving Related Settings

### Matrix Storage Scheme

#### Asymmetric Banded

By default, an **asymmetric banded** matrix storage scheme is used. For 1D analysis, the global stiffness matrix is\
always symmetric. However, when it comes to 2D and 3D analyses, the global stiffness may be structurally symmetric in\
terms of its layout, but there may not have the same value on the symmetric entries. In maths language, if $$K(i,j) \neq0$$ then $$K(j,i)\neq0$$, but $$K(i,j)\neq{}K(j,i)$$. Hence, an asymmetric banded storage is the safest.\
The `_gbsv()` LAPACK subroutine is used for matrix solving.

#### Symmetric Banded

It shall be noted that the symmetric scheme can save almost $$50\%$$ of the memory used in the asymmetric scheme.\
The `_pbsv()` LAPACK subroutine is used. This subroutine can only deal with **symmetric positive definite band matrix**.\
For some problems, in which the matrix is not necessarily positive definite, for example the buckling problems, this\
subroutine fails. **Before enabling the symmetric banded storage, the analyst must check if the matrix is SPD.**

#### Full Storage

If the problem scale is small, it does not hurt if a full storage scheme is used. For some particular issues such as\
particle collision issues, the full storage scheme is the only option.

#### Full Packed Storage

If the matrix is symmetric, a so-called pack format can be used to store the matrix. Essentially, only the upper or the\
lower triangle of the matrix is stored. The spatial cost is half of that of the full storage, but the solving speed is\
no better. The `_spsv()` subroutine is used for matrix solving. It is not recommended using this packed scheme.

#### Sparse Storage

The sparse matrix is also supported. Several sparse solvers are implemented.

### Direct System Solver

Different solvers are implemented for different storage schemes. It is possible to switch from one to another by using\
the following command. Details are covered in the summary table.

```
set system_solver (1)
# (1) string, system solver name
```

#### Mixed Precision Algorithm

The following command can be used to control if to use mixed precision refinement.\
This command has no effect if the target matrix storage scheme has no mixed precision implementation.

```
set precision (1)
# (1) string, "single" ("mixed") or "double" ("full")
```

#### Iterative Refinement

The mixed precision algorithm requires iterative refinement.\
The maximum number of refinements can be bounded by

```
set iterative_refinement (1)
# (1) integer, maximum number of refinements
```

It **cannot** exceed 256.

#### Iterative Refinement Tolerance

If the mixed precision algorithm is used, it is possible to use the following command to control the tolerance.

```
set tolerance (1)
# (1) double, tolerance
```

Typically, each refinement reduces the error by a factor of $$10^{-7}$$.\
Thus, two or three refinements should be sufficient to achieve the working precision.

Thus, the following command set makes sense.

```
set precision mixed
set iterative_refinement 3
set tolerance 1e-15
```

### Iterative System Solver

Iterative solvers are available for sparse storage.\
The [`Lis`](https://www.ssisc.org/lis/) library provides a wide range of solvers and preconditioners.\
See [here](https://github.com/TLCFEM/suanPan-manual/blob/dev/docs/Collection/Process/lis.md) for more details.

It is also possible to use GPU-based iterative solver powered by the [MAGMA](https://icl.utk.edu/magma/) library.\
The binaries shipped officially are **not** compiled with GPU support.\
Users can compile the library with GPU support by themselves.\
See [here](https://github.com/TLCFEM/suanPan-manual/blob/dev/docs/Collection/Process/magma.md) for more details.

### Summary

#### For Single Node Machine

With `SP_ENABLE_MPI` **disabled**, all available settings are summarised in the following table.

| storage       | configuration         | configuration        | system solver  | mixed precision | subroutine in external library |
| ------------- | --------------------- | -------------------- | -------------- | --------------- | ------------------------------ |
| full          | `set symm_mat false`  | `set band_mat false` | `LAPACK`       | yes             | (default) `d(s)gesv`           |
|               |                       |                      | `CUDA`         | yes             | `cusolverDnD(S)gesv`           |
| symm. banded  | `set symm_mat true`   | `set band_mat true`  | `LAPACK`       | yes             | (default) `d(s)pbsv`           |
|               |                       |                      | `SPIKE`        | yes             | `d(s)spike_gbsv`               |
| asymm. banded | `set symm_mat false`  | `set band_mat true`  | (not required) | yes             | `d(s)gbsv`                     |
| symm. packed  | `set symm_mat true`   | `set band_mat false` | (not required) | yes             | `d(s)ppsv`                     |
| sparse        | `set sparse_mat true` |                      | `SuperLU`      | no              | (default) `d(s)gssv`           |
|               |                       |                      | `CUDA`         | yes             | `cusolverSpD(S)csrlsvqr`       |
|               |                       |                      | `PARDISO`      | no              | `pardiso`                      |
|               |                       |                      | `FGMRES`       | no              | `dfgmres`                      |

#### For Multi Node Cluster

With `SP_ENABLE_MPI` **enabled**, all available settings are summarised in the following table.

| storage       | configuration         | configuration        | system solver  | mixed precision | subroutine in external library |
| ------------- | --------------------- | -------------------- | -------------- | --------------- | ------------------------------ |
| full          | `set symm_mat false`  | `set band_mat false` | (not required) | no              | `pdgesv`                       |
| symm. banded  | `set symm_mat true`   | `set band_mat true`  | (not required) | no              | `pdpbsv`                       |
| asymm. banded | `set symm_mat false`  | `set band_mat true`  | (not required) | no              | `pdgbsv`                       |
| symm. packed  | `set symm_mat true`   | `set band_mat false` | (not required) | no              | `pdposv`                       |
| sparse        | `set sparse_mat true` |                      | `PARDISO`      | no              | `cluster_sparse_solver`        |
|               |                       |                      | `LIS`          | no              | `lis_solve`                    |
|               |                       |                      | `MUMPS`        | no              | (default) `dmumps_c`           |

Some empirical guidance can be concluded as follows.

1. For most cases, the asymmetric banded storage with full precision solver is the most general option.
2. The best performance is obtained by using symmetric banded storage, if the (effective) stiffness matrix is guaranteed\
   to be positive definite, users shall use it as a priority.
3. The mixed precision algorithm often gives the most significant performance boost for full storage with `CUDA` solver.\
   It outperforms the full precision algorithm when the size of system exceeds several thousands.
4. The `SPIKE` solver is slightly slower than the conventional `LAPACK` implementations.
5. The `PARDISO` direct solver and `FGMRES` iterative solver are provided by `MKL`.

## Parallel Matrix Assembling

For dense matrix storage schemes, the global matrix is stored in a consecutive chunk of memory. Assembling global matrix\
needs to fill in the corresponding memory location with potentially several values contributed by different elements.\
Often in-place atomic summation is involved. To assign each memory location a mutex lock is not cost-efficient. Instead,\
a $$k$$-coloring concept can be adopted to divide the whole model into several groups. Each group contains elements that\
do not share common nodes with others in the same group. By such, no atomic operations are required. Elements in the\
same group can be updated simultaneously so the matrix assembling is lock free.

By default, the coloring algorithm is enabled. To disable it, users can use the following command.

```
set color_model false
set color_model none
```

As a NP hard problem, there is no optimal algorithm to find the minimum chromatic number. The Welsh-Powell algorithm is\
implemented in `suanPan`. The maximum independent set algorithm is also available, it may outperform the\
Welsh-Powell algorithm on large models. To switch, users can use the following command.

```
# default to WP algorithm
set color_model WP
# switch to MIS algorithm
set color_model MIS
```

Also, depending on the problem setup, such a coloring may or may not help to improve the performance. If there are not a\
large number of matrix assembling, the time saved may not be significant. Thus, for problems of small sizes, users may\
consider disabling coloring.

This option has no effect if a sparse storage is used.

## Penalty Number

For some constraints and loads that are implemented by using the penalty method, the default penalty number can be\
overridden.

```
set constraint_multiplier (1)
set load_multiplier (1)
# (1) double, new penalty number
```

This command does not overwrite user defined penalty number if the specific constraint or load takes the penalty number\
than input arguments.

## FGMRES Iterative Tolerance

For `FGMRES` iterative solver, one can use the following dedicated command to control the tolerance of the algorithm.

```
set fgmres_tolerance (1)
# (1) double, tolerance
```

If the boundary condition is applied via the penalty method, say for example one previously\
uses `set constraint_multiplier 1E8`, then there is no need to set a tolerance smaller than `1E-8`. A slightly larger\
value is sufficient for an iterative algorithm, one can then set

```
set fgmres_tolerance 1E-6
```
