Tutorial 5b – Inhomogeneous one-dimensional log-Normal random field

#! ===============================
#! Load the FE module
#! ===============================

loadlib "fem";

#! ===============================
#! Initial definitions
#! ===============================

const l  = 3;              # [m] length of random field domain
const mu = 10;             # mean value of random field
var sd = 4*(gx-1.5)^2+1;   # std. dev. function of random field

# auto-correlation coefficient function of underlying Gaussian field
var corr = exp(-1*(deltap/(min((min(gx,gx2)/1.5)^4,2)+0.1))^2);

const M  = 10;     # number of terms in the random field expansion

# mean of underlying Gaussian field
  fun lambda(2) = ln($1)-ln(($2/$1)^2+1)/2;
# std. dev. of underlying Gaussian field
  fun chi(2)    = sqrt(ln(($2/$1)^2+1));


#! ===============================
#! Define the mesh
#! ===============================
# The mesh becomes finer from right to left. 

const Nel = 20;
const x = 1;
fun xf(2) = $1*(1-1/$2)^(($2/Nel)^8*1+1);
for (i=Nel;i>0.5;i-1) {
  node i [l*x];
  funplot i, l*x;
  const x = xf(x,i);
};
node 0 [0];
for (i=1;i<Nel+0.5;i+1) {
  edge i i-1 i;
};


#! ===============================
#! Random field discretization (FE-KL)
#! ===============================

const pO = 10;   # polynomial order of the FE shape functions
const MO = 0.1;  # 1/MO is number of sub-intervals in outer integration loop

# type of discretization: finite element KL expansion
  rf_new 1 kl_fem_gauss mu lambda(mu,sd) sd chi(mu,sd) rhocor CORR;
# Assign a geometry to the field:
  for (i=1;i<Nel+0.5;i+1) {
    rf_elset 1 edge i porder=PO MAXSIZEO=MO {dolog=false;}; 
  };
# Obtain the random field discretization
  rf_solve 1 dim=M calcerr=false; 

funplot #rf_err(1,1), 
  rf_err(1,101),      # Var-based measure - used instead of '1' as std.dev. varies
  rf_err(1,2),        # Cov-based measure
  rf_err(1,3);        # error measure 3

#! ===============================
#! Output eigenvalues
#! ===============================

const NEV = M;
sfor (i;NEV) {
  funplot i, rf_eigenvalue(1,i);
};

#! ===============================
#! Plot the eigenmodes of the random fields
#! ===============================

#! -------------------------------
#! output the eigenmodes to files
#! -------------------------------

const ppd = 500;                # number of points to use for plotting
mtxconst_new ppd1(3,1,ppd/Nel); # ... for 2nd random field

filestream fs ("output_rf2.dat");
rf_plot 1 { ppdim = ppd1; stream=fs; };


#! -------------------------------
#! invoke the tool matplotlib
#! -------------------------------
# You can comment this section out if you do not have 
# 'Python' or 'matplotlib' installed on your system.

# load the python-module
  loadlib "Python";
# add current path to Python search path
  python_path $pwd();
# load the module 'tutorial_1d_lsf'
  python_module tutorial_5b_plot;
# link the Python function to Fesslix
  python_fun plot_fun = tutorial_5b_plot.PlotModes; 
# start the plotting 
  calc plot_fun();

#! ===============================
#! Plot covariance around 1.5
#! ===============================
# of the underlying Gaussian random field!!!

#! -------------------------------
#! output the covariance structure
#! -------------------------------

filestream fs ("output_rf_cov.dat");
rf_plotcov 1 [1.5] {stream=fs;};
  
#! -------------------------------
#! invoke the tool matplotlib
#! -------------------------------
# You can comment this section out if you do not have 
# 'Python' or 'matplotlib' installed on your system.

# link the Python function to Fesslix
  python_fun plot_cov = tutorial_5b_plot.PlotCov; 

calc plot_cov();

#! ===============================
#! Plot realizations
#! ===============================

#! -------------------------------
#! generate realizations
#! -------------------------------
const Nrep = 10;
const Npel = 20;
mtxconst_new resv(Nel*(Npel+1),Nrep+1,0);
sfor (j;Nrep) {
  const c = 0;
  rnd_smp "rf_1";
  for (i=1;i<Nel+0.5;i+1) {
    for (lx=-1;lx<=1;lx+2/Npel) {
      const res = geocalc(edge,i,exp(rf(1)),lx);
      mtxcoeff resv(c,0) = gx;
      mtxcoeff resv(c,j) = res;
      const c += 1;
    };
  };
};

#! -------------------------------
#! output the realizations
#! -------------------------------
filestream fs("output_reali.dat");
sfor (j;Nel*(Npel+1)) {
  mtxconst_sub vsub = resv(row:j);
  funplot {vsub} { stream=fs; };
};

#! -------------------------------
#! invoke the tool matplotlib
#! -------------------------------
# You can comment this section out if you do not have 
# 'Python' or 'matplotlib' installed on your system.

# link the Python function to Fesslix
  python_fun plot_reali = tutorial_5b_plot.PlotReali; 

calc plot_reali();