statROM/demo_2d.py at master · herluc/statROM · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
"""
statROM Helmholtz 2D scattering example

This is the high level run-file. User parameters are declared and the individual parts
of the statROM procedure are called. exampleHelmholtz.py is imported to provide the relevant
methods for data generation, prior computation, error estimation and data assimilation.
"""

__author__ = "Lucas Hermann"
__version__ = "0.1.1"
__license__ = "MIT"


import numpy as np
from exampleHelmholtz_scatter import StatROM_2D


class UserParameters:
    def __init__(self):
        pass
    f = 400#Hz      # frequency
    s = 250#Hz      # AORA expansion frequency
    m = 12#8#12          # number of matched moments in AORA

    #               # the mesh for the model and a finer one for data are hardcoded with gmsh.
    #               #
    n_qmc = 256     # number of QMC sample points. Only powers of two are permitted

    ns = 5          # number of sensors
    no = 20         # number of observations per sensor
    sig_o = 5e-4#Pa # observation noise

    n_est = 200     # number of error estimator training points


if __name__ == '__main__':
    up = UserParameters
    funcs = StatROM_2D(up)

    ### Generate data from fine mesh:
    funcs.switchMesh_self("ground_truth")
    funcs.generateParameterSamples(up.n_qmc,save=False,load=False, seed=8) # seed = None for a new sample at each new run
    #funcs.generateParameterSamples(up.n_qmc,save=True,load=False)
    funcs.getFullOrderPrior() # solve for "ground truth"
    funcs.get_noisy_data_from_solution(0,np.abs(funcs.u_mean_std)) # extract noisy data at sensors
    ####

    funcs.switchMesh_self("coarse")

    funcs.generateParameterSamples(up.n_qmc,load = False, seed=9) # get new parameter sample, seed = None for a new sample at each new run

    # Full order reference:
    funcs.getFullOrderPrior()

    ### Compute ROM basis:
    _, funcs.V_adj_list = funcs.getAORAbasis(Nr=funcs.L,rhs_sp=np.zeros(np.shape(funcs.RBmodel.coordinates_coarse)[0]),adj=True)[0:2] # adjoint solves
    for i,sample in enumerate(funcs.f_samples):
        print("primal sample number: "+str(i))
        funcs.computeROMbasisSample(sample,i)
    ####

    # ROM prior and error estiamte:
    funcs.calcROMprior()

    funcs.plotPriorVtk()

    # FOM reference posterior, classical ROM posterior and proposed ROM posterior:
    funcs.getFullOrderPosterior()
    funcs.getEasyROMPosterior()
    funcs.getCorrectedROMPosterior()

    funcs.plotPosteriorVtk()