Tutorial cpt random fields

docs/tutorials.rst

@@ -984,4 +984,308 @@ The calculation is then ran by calling the run_calculation function within the s
 
 .. code-block:: python
 
-    stem.run_calculation()
+    stem.run_calculation()
+
+
+.. _tutorial4:
+
+CPT-based random field model with point load
+--------------------------------------------
+This tutorial shows how a subsurface model can be created directly from CPT data. First, the concept of the CPT-based random field generation for subsurface models is presented. Secon
+
+
+Theory
+......
+
+* CPT data is read and interpreted into a set of point with coordinates `(x,y,z)` and corresponding interpreted properties $z_{data}$ (e.g. shear wave velocity `vs` based on  tip resistance and sleeve friction). The the distribution of the dataset $x_{data}$ is considered to be representative for the distribution of the site $x\inX$. 
+* A transfomration model is created based on the CPT data. This transfomration model $X=T(Z)$ characterises the distribution of the physical stochastic parameter $X$ as a function of the standard-normal stochastic variable $Z \\sim N(0,1)$. This model is used to create the standard-normal equivalent data $z_{data} = T^{-1}(x_{data})$ 
+* The standard-normal data is used to calibrate a Gaussian regression model. Sampling from this model provides the conditioned . Internally, the sampoling is performed in the traditional Kriging formulation, in which an unconditioned random field is conditioned to the conditioning data . This takes place in the standard-normal space, resultin gin standar-normal fields $u_{field}$.
+* The generated conditioned random fields are transformed to the parameter space using the transformation model: $v_{field} = T(u_{field})$.
+
+
+Implementation
+..............
+The CPT-based generation is implemented in a module of the `stem.random_fields` package. A single wrapper class is provided to combine the different components into a single random field generator that can be linked to STEM.
+
+
+.. code_block:: python
+
+    from geostatistical_cpt_interpretation import CPT_data
+
+    cpt_folder = r'/benchmark_tests/test_cpt_conditioning/cpt_data'
+    cpt_data = CPT_data(cpt_directory = cpt_folder)
+    cpt_data.read_cpt_data()
+    cpt_data.interpret_cpt_data()
+
+
+    # visualise
+    plt.title('Original coordinates')
+    cpt_data.plot_coordinates()
+
+    # transform the data to a modelling domain
+    cpt_data.data_coordinate_change(orientation_x_axis=73,based_on_midpoint=True)
+
+    plt.title('Model coordinates')
+    cpt_data.plot_coordinates()
+
+
+The `CPT_data` class makes use of the `d-geolib-plus` package, and contains all CPT data after reading. If necessary, the data is accessible, for example: 
+
+.. code-block:: python
+
+    from pathlib import Path
+    # create plots of CPT data using geolib-plus functionalities 
+    for cpt in cpt_data.cpt_list:
+        cpt.plot(Path('.'))
+
+The marginal transformation model can be generated and visualised:
+
+.. code-block:: python
+
+    marginal_transformator = MarginalTransformation(cpt_data.vs,min = 50)
+    marginal_transformator.plot(x_label = '$u$ : standard-normal variable',y_label = '$v$ : shear wave velocity [m/s]')
+
+
+Next, the data for the calibration of the geostatistical model can be selected. To allow a faster calibration of the geostatistical model, only a selection of the data is used. This selection of 2000 pooints is made randomly. The data is transformed to standard-normal data using the transformator and only the `z` (horizontal) and `y` (vertical) coordinates are selected. Calibration is based on likelihood maximisation.   
+
+.. code-block:: python
+
+    index_selection = np.random.choice(len(cpt_data.vs),size = 2000,replace = False)
+
+    coords = cpt_data.data_coords[index_selection]
+    z_data = marginal_transformator.x_to_z(x = cpt_data.vs[index_selection])
+
+    geo_model = GeostatisticalModel(nb_dimensions=2,v_dim = 1)
+    geo_model.calibrate(coords = coords[:,[2,1]],values = z_data)
+
+Next, the random field properties are transferred to the random field generator. This generator needs to be initiated with the same spatial correlation model as used for the calibration (default = Gaussian). In addition, conditioning points are required to generate meaningful conditioned random fields. These points can, but don't need to be, the same points as used for the calibration. Also, conditional simulation can account for the noise in the calibration and conditioning data. Thios noise is not included in the field itself, but instead allows for a small deviation of the generated random fields from the conditioning points:
+
+.. code-block:: python
+
+    random_field_generator = RandomFields(model_name = ModelName.Gaussian, 
+                                n_dim = 2, 
+                                mean = 0, 
+                                variance = 1,
+                                v_scale_fluctuation = geo_model.vertical_scale_fluctuation, 
+                                anisotropy = geo_model.anisotropy, 
+                                angle = [0], 
+                                seed = 13)
+
+    I = np.random.choice(len(cpt_data.vs),size = 3000,replace = False)
+
+    coords = cpt_data.data_coords[I]
+    values = cpt_data.vs[I]
+    random_field_generator.set_conditioning_points(points = coords[:,[2,1]],
+                                values = marginal_transformator.x_to_z(x = values),
+                                noise_level = geo_model.noise_level)
+
+A conditioned random field is generated on a regular grid of coordinates:
+
+.. code-block:: python
+
+    # create grit of points on the domnain (-220,220) by (-24,-1) to generate a field for.
+    x = np.linspace(-220,220,250)
+    z = np.linspace(-24,-1,250)
+    X,Z = np.meshgrid(x,z)
+
+    # generate a conditioned random field
+    sample_coords = np.array([X.ravel(),Z.ravel()]).T
+    random_field_generator.generate_conditioned(nodes = sample_coords)
+
+The conditioned random field that is generated contains values at the prediction points (the regular grid) as well as on the conditioning point coordinates. In the generated array, the values generated at the conditioning point coordinates are at the end. They are generally not needed as part of the generated random field and can be left out. The generated standard-normal field needs to be transformed to the marginal distribution of the physical variable (in this case the shear wave velocity) by the marginal transformation   
+
+.. code-block:: python
+
+    # Transform the generated standard-normal field to the distribution of the shear wave velocity 
+    z_map = random_field_generator.conditioned_random_field
+    vs_map = marginal_transformator.z_to_x(z_map[:250*250].reshape([250,250]))
+
+    # visualisation usimg `matplotlib.pyplot as plt`
+    plt.contourf(X,Z,vs_map)
+    plt.scatter(cpt_data.data_coords[:,2],cpt_data.data_coords[:,1],c = cpt_data.vs)
+
+Finally, visualisation gives the following 2D cross-section along 400 meters of the Delft-Schiedam line
+
+
+
+.. image:: _static/2D_conditional_random_field_vs.png
+
+
+[INTRODUCE DIFFERENT CLASSES HERE]
+
+
+
+Application in STEM
+...................
+
+The modules required for the computational model are loaded, together with two additional classes: `stem.additional_processes.ParamaterFieldParameter` is used as the interface between the `STEM` package and the `random_fields` package; `random_fields.geostatistical_cpt_interpretation.ElasticityFieldsCpt` is used as a generator of conditioned random fields.
+
+.. code-block:: python
+
+    import os
+    from stem.model import Model
+    from stem.soil_material import OnePhaseSoil, LinearElasticSoil, SoilMaterial, SaturatedBelowPhreaticLevelLaw
+    from stem.load import PointLoad
+    from stem.boundary import DisplacementConstraint
+    from stem.solver import AnalysisType, SolutionType, TimeIntegration, DisplacementConvergenceCriteria,\
+         NewtonRaphsonStrategy, NewmarkScheme, Amgcl, StressInitialisationType, SolverSettings, Problem
+    from stem.output import NodalOutput, VtkOutputParameters, Output, GaussPoi    ntOutput
+    from stem.stem import Stem
+
+    from stem.additional_processes import ParameterFieldParameters
+    from geostatistical_cpt_interpretation import ElasticityFieldsFromCpt
+
+As an example, a model for a single block of soil is created:
+
+.. code-block:: python
+
+    ndim = 3
+    model = Model(ndim)
+
+    solid_density_1 = 2650
+    porosity_1 = 0.3
+    young_modulus_1 = 0.
+    poisson_ratio_1 = 0.2
+    soil_formulation_1 = OnePhaseSoil(ndim, IS_DRAINED=True, DENSITY_SOLID=solid_density_1, POROSITY=porosity_1)
+    constitutive_law_1 = LinearElasticSoil(YOUNG_MODULUS=young_modulus_1, POISSON_RATIO=poisson_ratio_1)
+    retention_parameters_1 = SaturatedBelowPhreaticLevelLaw()
+    material_soil_1 = SoilMaterial("soil_1", soil_formulation_1, constitutive_law_1, retention_parameters_1)
+
+    soil1_coordinates = [( 0.0, -25.0, -25.0), 
+                         ( 20.0, -25.0, -25.0), 
+                         ( 20.0,  -1.0, -25.0), 
+                         ( 1.0,  -1.0, -25.0), 
+                         ( 0.0,  -1.0, -25.0)]
+
+    model.extrusion_length = 50.
+    orientation_x_axis = 72.
+
+    model.set_mesh_size(element_size=1.)
+
+    model.add_soil_layer_by_coordinates(soil1_coordinates, material_soil_1, "soil_layer_1")
+
+    # create input files directory, since it might not have been created yet
+    os.makedirs(input_files_dir, exist_ok=True)
+
+The random field generator for the Young modulus is set up as follows, using only `y` and 'z' coordinates for the calibration, because the CPTs are all on one line along the z-axis. The length scale for the z-axis will be used for the x-axis as well. All other settings are identical to earlier examples:
+
+.. code-block:: python
+
+    cpt_folder = r'/benchmark_tests/test_cpt_conditioning/cpt_data'
+    orientation_x_axis = 75
+
+    elastic_field_generator_cpt = ElasticityFieldsFromCpt(cpt_file_folder = cpt_folder,
+                based_on_midpoint = True,
+                max_conditioning_points = 1000,
+                orientation_x_axis = orientation_x_axis,
+                return_property = 'young_modulus')
+    elastic_field_generator_cpt.calibrate_geostat_model(calibration_indices=(1,2),v_dim=0)
+
+    field_parameters_json = ParameterFieldParameters(
+                property_name="YOUNG_MODULUS",
+                function_type="json_file",
+                field_generator=elastic_field_generator_cpt)
+
+    # add the random field to the model
+    model.add_field(part_name="soil_layer_1", field_parameters=field_parameters_json)
+
+Add a line load at the boundary. The mesh around the point load is [TODO: FIX ELEMENT SIZE AT LINE LOAD: set_element IS NOT ACCOUNTED FOR IN THE CODE]:
+
+.. code-block:: python
+
+    load_coordinates = [(1., -1.0, 25),(1., -1.0, 25)]
+    point_load = PointLoad(active=[False, True, False], value=[0, -10000, 0])
+    model.add_load_by_coordinates(load_coordinates, point_load, "line_load")
+
+    model.set_element_size_of_group(0.5,'line_load')
+
+    model.show_geometry(show_surface_ids=True)
+
+Add boundary conditions [TODO: FIX TO REALISTIC BOUNDARY CONDITIONS]:
+
+.. code-block:: python
+
+    no_displacement_parameters = DisplacementConstraint(active=[True, True, True],
+                                                        is_fixed=[True, True, True], value=[0, 0, 0])
+    roller_displacement_parameters = DisplacementConstraint(active=[True, True, True],
+                                                            is_fixed=[True, False, True], value=[0, 0, 0])
+
+
+    model.add_boundary_condition_by_geometry_ids(2, [2], no_displacement_parameters, "base_fixed")
+    model.add_boundary_condition_by_geometry_ids(2, [1, 7], roller_displacement_parameters, "sides_roller")
+
+
+Set problem [TODO: CHANGE TO QUASI STATIC IN A SINGLE STEP]: 
+
+.. code-block:: python
+
+    end_time = 1.
+    delta_time = 1.
+    analysis_type = AnalysisType.MECHANICAL
+    solution_type = SolutionType.DYNAMIC
+
+    # Set up start and end time of calculation, time step and etc
+    time_integration = TimeIntegration(start_time=0.0, end_time=1., delta_time=0.25, reduction_factor=1.0,
+                                       increase_factor=1.0)
+    convergence_criterion = DisplacementConvergenceCriteria(displacement_relative_tolerance=1.0e-4,
+                                                            displacement_absolute_tolerance=1.0e-9)
+    strategy_type = NewtonRaphsonStrategy()
+    scheme_type = NewmarkScheme()
+    linear_solver_settings = Amgcl()
+    stress_initialisation_type = StressInitialisationType.NONE
+    solver_settings = SolverSettings(analysis_type=analysis_type, solution_type=solution_type,
+                                     stress_initialisation_type=stress_initialisation_type,
+                                     time_integration=time_integration,
+                                     is_stiffness_matrix_constant=True, are_mass_and_damping_constant=True,
+                                     convergence_criteria=convergence_criterion,
+                                     strategy_type=strategy_type, scheme=scheme_type,
+                                     linear_solver_settings=linear_solver_settings, rayleigh_k=0.12,
+                                     rayleigh_m=0.0001)
+
+
+    problem = Problem(problem_name="calculate_load_on_spatially_variable_embankment_3d", number_of_threads=1,
+                      settings=solver_settings)
+    model.project_parameters = problem
+
+
+Define output. Note that to visualise the random fields, here the Young modulus is set as a variable to include in the output. This is done at the integration points:
+
+.. code-block:: python
+    nodal_results = [NodalOutput.DISPLACEMENT, 
+                     NodalOutput.VELOCITY, 
+                     NodalOutput.ACCELERATION]
+    gauss_point_results = [GaussPointOutput.YOUNG_MODULUS]
+
+    model.add_output_settings(
+        part_name="porous_computational_model_part",
+        output_dir=results_dir,
+        output_name="vtk_output",
+        output_parameters=VtkOutputParameters(
+            file_format="ascii",
+            output_interval=1,
+            nodal_results=nodal_results,
+            gauss_point_results=gauss_point_results,
+            output_control_type="step"
+        )
+    )
+
+The code is run in a single push as before:
+
+.. code-block:: python
+
+    stem = Stem(model, input_files_dir)
+    stem.write_all_input_files()
+    stem.run_calculation()
+
+
+Visualisation in Paraview provides the following view of Young modulus values at 3 of the 5 CPTs and a block of soil in which the conditioned random field represents a field of Young modulus values in the domain:
+
+.. image:: _static/3D_conditioned_random_field_young_modulus.PNG
+
+
+
+
+
+
+
+