Module reducedbasis

1. Import

The work has not yet been merged in develop. To include the module at current state, you need to checkout the branch feature/pyrb.

Line to import the module :

import feelpp.mor.reducedbasis.reducedbasis_timeOffline as mor_rb
This import load the module for offline generation, associated to the time-dependent case. For the time-independent case, you can import feelpp.mor.reducedbasis.reducedbasisOffline instead.

To set the environment, those module also need to be imported :

import sys, os
from feelpp.toolboxes.heat import * (1)
import feelpp.toolboxes.core as core (2)
from feelpp.operators import mass (3)
import feelpp.mor as mor (4)
import feelpp
1 The toolbox heat is used to simulate the heat equation on the case studied.
2 The module core is used to set the environment.
3 The mass function is needed to create the mass matrix if we work on a time-dependent case.
4 The module mor is used to get the functions neededd by the toolbox mor.

2. Set the environment

config = feelpp.localRepository("reducedbasis")
sys.argv = ["reducedbasis"]
o = core.toolboxes_options("heat")
o.add(mor.makeToolboxMorOptions())

e = feelpp.Environment(sys.argv, opts=o, config=config)

3. Set the model

3.1. Description of the model

In this document, the cases used is thermal-fin in 2D, with a \(\mathbb{P}_1\) discretization and stationnary. This model is described in this page.

The usage would be the same for another case.

3.2. Code

Set the parameters nedded to the model :

  • config_file : path to the cfg file

  • dim : the dimension of the case

  • order : the order of discretization (1 for \(\mathbb{P}_1\), 2 for \(\mathbb{P}_2\))

  • time_depedent : is the case stationnary or transient ? (for now always to False)

  • name : name of the CRB case

config_file = "thermal-fin.cfg"
dim = 2
order = 1
time_depedent = False
name = "thermalfin-2d"

The configuration files for this case are available on github : geometry file, cfg file, and json file.

Then, we can create the model :

cfg = feelpp.readCfg(config_file)
feelpp.Environment.setConfigFile(config_file)
model_path = feelpp.Environment.expand(cfg['heat']['filename'])
j = feelpp.readJson(model_path)
try:
    j.pop('PostProcess')
except KeyError as e:
    print(f"There was no section {e} in the model")

Now we load the crb parameters, the path has to be set as the option toolboxmor.filename in the cfg file.

crb_model_file = feelpp.Environment.expand(cfg['toolboxmor']['filename'])
crb_model_file
crb_model_properties = mor.CRBModelProperties(worldComm=feelpp.Environment.worldCommPtr())
crb_model_properties.setup(crb_model_file)
crb_model_outputs = crb_model_properties.outputs()

output_names = []
for n, _ in crb_model_outputs:
    output_names.append(n)

print(f"Outputs of the models are {output_names}")

4. Set the toolboxes

heatBox = heat(dim=dim, order=order)
heatBox.init()

The parameter default_parameter is defined in the section Parameters of the JSON file.

modelParameters = heatBox.modelProperties().parameters()
default_parameter = modelParameters.toParameterValues()
model = mor.toolboxmor(name=name, dim=dim, time_dependent=time_depedent)
model.setFunctionSpaces( Vh=heatBox.spaceTemperature() )

def assembleDEIM(mu):
    for i in range(0,mu.size()):
        heatBox.addParameterInModelProperties(mu.parameterName(i), mu(i))
    heatBox.updateParameterValues()
    return heatBox.assembleRhs()

def assembleMDEIM(mu):
    for i in range(0,mu.size()):
        heatBox.addParameterInModelProperties(mu.parameterName(i), mu(i))
    heatBox.updateParameterValues()
    return heatBox.assembleMatrix()

model.setAssembleDEIM(fct=assembleDEIM)
model.setAssembleMDEIM(fct=assembleMDEIM)

model.initModel()

Set the toolboxes associated to the reduced meshes :

heatBoxDEIM = heat(dim=dim, order=order)
heatBoxDEIM.setModelProperties(j)
meshDEIM = model.getDEIMReducedMesh()
heatBoxDEIM.setMesh(meshDEIM)
heatBoxDEIM.init()

heatBoxMDEIM = heat(dim=dim, order=order)
heatBoxMDEIM.setModelProperties(j)
meshMDEIM = model.getMDEIMReducedMesh()
heatBoxMDEIM.setMesh(meshMDEIM)
heatBoxMDEIM.init()

def assembleOnlineDEIM(mu):
    for i in range(0,mu.size()):
        heatBoxDEIM.addParameterInModelProperties(mu.parameterName(i),mu(i))
    heatBoxDEIM.updateParameterValues()
    return heatBoxDEIM.assembleRhs()

def assembleOnlineMDEIM(mu):
    for i in range(0,mu.size()):
        heatBoxMDEIM.addParameterInModelProperties(mu.parameterName(i),mu(i))
    heatBoxMDEIM.updateParameterValues()
    return heatBoxMDEIM.assembleMatrix()

model.setOnlineAssembleDEIM(assembleOnlineDEIM)
model.setOnlineAssembleMDEIM(assembleOnlineMDEIM)

model.postInitModel()
model.setInitialized(True)

5. Offline generation

Information about the parameter space \(D^\mu\)

Dmu = model.parameterSpace()
print(Dmu.parameterNames())
print(Dmu.mumin())
print(Dmu.mumax())
Results 
['Bi'    ,'k_1'   ,'k_2'   ,'k_3'   ,'k_4']
[1.00e-02,1.00e-01,1.00e-01,1.00e-01,1.00e-01]
[1.00e+00,1.00e+01,1.00e+01,1.00e+01,1.00e+01]
mubar = Dmu.element(True, False)
mubar.setParameters(default_parameter)
print("mubar =")
mubar.view()
Results 
mubar =
Bi	0.01
k_1	0.1
k_2	0.1
k_3	0.1
k_4	0.1

Get affine decomposition \(A(\mu) = \displaystyle\sum_q \beta^q_A(\mu) A^q\), \(F(\mu)=\displaystyle\sum_p \beta_F^p(\mu) F^(\mu)\). To use the matrices and vector in the class, we need to convert themn to petsc4py objects, using the functions convertToPetscMat and convertToPetscVec.

affineDecomposition = model.getAffineDecomposition()
Aq_ = affineDecomposition[0]
Fq_ = affineDecomposition[1]

Aq = mor_rb.convertToPetscMat(Aq_[0])
Fq = []
for f in Fq_:
    Fq.append(mor_rb.convertToPetscVec(f[0]))

print("Aq", Aq)
print("Fq", Fq)
Results 
Aq [<petsc4py.PETSc.Mat object at 0x7f5f693abba0>, <petsc4py.PETSc.Mat object at 0x7f5f7000d210>, <petsc4py.PETSc.Mat object at 0x7f5f693fd3f0>, <petsc4py.PETSc.Mat object at 0x7f5f693fd3a0>, <petsc4py.PETSc.Mat object at 0x7f5f693fd210>, <petsc4py.PETSc.Mat object at 0x7f5f693fd260>]
Fq [[<petsc4py.PETSc.Vec object at 0x7f5f693fd350>], [<petsc4py.PETSc.Vec object at 0x7f5f693fd2b0>], [<petsc4py.PETSc.Vec object at 0x7f5f693fd1c0>]]

Now we create the object reduced_basis.

rb = mor_rb.reducedbasisOffline(Aq=Aq, Fq=Fq, model=model, mubar=mubar,
                output_names=output_names, use_dual_norm=False)
Results 
[reducedbasis] Online rb initialized
[slepc4py] number of (smaller) eigenvalues computed : 16
[reducedbasis] Constant of coercivity : 0.9999999999999977
[slepc4py] number of eigenvalues computed : 16
[reducedbasis] Constant of continuity : 1.0000000000000018
[reducedbasis] Offline rb initialized

Generation of the basis, using the greedy algorithm.

Ntrain = 1000
s = Dmu.sampling()
s.sampling(Ntrain, "random")
Xi_train = s.getVector()
mu_0 = Dmu.element()
rb.greedy(mu_0, Xi_train, eps_tol=1e-6)
Results 
[reducedBasis] Start greedy algorithm
[reducedBasis] Computing betas: 100%|█████████████████████████████████████████████████| 999/999 [00:25<00:00, 39.91it/s]
[reducedBasis] Offline generation of the basis: 100%|█████████████████████████████████████| 1/1 [00:00<00:00, 32.94it/s]
[reducedBasis] Gram-Schmidt orthonormalization done after 1 step
[reducedBasis] Greedy, step 1: 100%|████████████████████████████████████████████████| 999/999 [00:00<00:00, 3301.45it/s]
[reducedBasis] Greedy algo, N=1, Δ=2.735059e-06 (tol=1.000000e-06) µ=[2.17e-01,8.67e+00,7.04e+00,4.51e+00,8.96e-01]
[reducedBasis] Gram-Schmidt orthonormalization done after 1 step
[reducedBasis] Greedy, step 2: 100%|████████████████████████████████████████████████| 998/998 [00:00<00:00, 3034.75it/s]
[reducedBasis] Greedy algo, N=2, Δ=5.947858e-07 (tol=1.000000e-06) µ=[4.39e-01,4.62e+00,4.23e-01,9.22e+00,4.94e-01]
[reducedBasis] End greedy algorithm

Compute offline errors :

rb.computeOfflineErrorRhs()
rb.computeOfflineError()

The script generate_basis.py run this offline part. See the dedicated page for more details.

6. Online phase

mu = Dmu.element()
mu.view()
Results 
Bi	0.851312
k_1	8.50443
k_2	0.687991
k_3	8.1523
k_4	3.47058

The function getSolutions returns the solutions of the model and the reduced model, the parameter k gives the index of the output wanted, -1 standing for the compliant output.

rb.getSolutions(mu, k=0)
Results 
(array([1.75329542e-03, 8.16719018e-10]), 3.2730508127626815e-06)

We can compare it to the finite element solution :

rb.getSolutionsFE(mu, k=0)
Results 
(<petsc4py.PETSc.Vec at 0x7f5f3d97eac0>, 3.273050812763351e-06)

The two functions getSolutions (resp. getSolutionsFE) return the solution \(u_N(\mu)\) (resp. \(u(\mu)\)) and the value of the output \(s_N^k(\mu)\) (res. \(s^k(\mu)\)). See the dedicated page for the API of ParameterSpaceElement.