Skip to content
Snippets Groups Projects
Verified Commit 99cd534d authored by Lars Bilke's avatar Lars Bilke
Browse files

[web] Added video tutorial by Dominik and Christian.

parent 2e0600f1
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 6250 additions and 0 deletions
web/content/docs/tools/getting-started/video-tutorial/OGSinput_basin.prj 0 → 100644
+ 212
− 0
View file @ 99cd534d
<?xml version="1.0" encoding="ISO-8859-1"?>
<OpenGeoSysProject>
<meshes>
<mesh>mesh_basin_domain.vtu</mesh>
<mesh>mesh_basin_physical_group_Left.vtu</mesh>
<mesh>mesh_basin_physical_group_Right.vtu</mesh>
<mesh>mesh_basin_physical_group_Bottom.vtu</mesh>
<mesh>mesh_basin_physical_group_Top.vtu</mesh>
</meshes>
<python_script>timeBCs_glacier.py</python_script>
<processes>
<process>
<name>SD</name>
<type>SMALL_DEFORMATION</type>
<!--define the numerical integration order -->
<!--(polynomial degree to which exact integration is possible) -->
<integration_order>2</integration_order>
<!--define the constitutive behavior with regard to the MaterialID -->
<constitutive_relation id="0">
<!-- soft sediment layer -->
<type>LinearElasticIsotropic</type>
<youngs_modulus>YoungModulus</youngs_modulus>
<poissons_ratio>PoissonRatio</poissons_ratio>
</constitutive_relation>
<constitutive_relation id="1">
<!-- stiff sediment layer -->
<type>LinearElasticIsotropic</type>
<youngs_modulus>YoungModulus</youngs_modulus>
<poissons_ratio>PoissonRatio</poissons_ratio>
</constitutive_relation>
<constitutive_relation id="2">
<!-- soft sediment layer -->
<type>LinearElasticIsotropic</type>
<youngs_modulus>YoungModulus</youngs_modulus>
<poissons_ratio>PoissonRatio</poissons_ratio>
</constitutive_relation>
<constitutive_relation id="3">
<!-- very stiff rock bed -->
<type>LinearElasticIsotropic</type>
<youngs_modulus>YoungModulus</youngs_modulus>
<poissons_ratio>PoissonRatio</poissons_ratio>
</constitutive_relation>
<solid_density>rho_sr</solid_density>
<specific_body_force>0 -9.81</specific_body_force>
<process_variables>
<process_variable>displacement</process_variable>
</process_variables>
<!--define the output quantities derived from the primary variables-->
<secondary_variables>
<secondary_variable internal_name="sigma" output_name="sigma"/>
<secondary_variable internal_name="epsilon" output_name="epsilon"/>
</secondary_variables>
</process>
</processes>
<time_loop>
<processes>
<process ref="SD">
<nonlinear_solver>basic_newton</nonlinear_solver>
<convergence_criterion>
<type>DeltaX</type>
<norm_type>INFINITY_N</norm_type>
<abstol>1e-14</abstol>
</convergence_criterion>
<time_discretization>
<type>BackwardEuler</type>
</time_discretization>
<time_stepping>
<type>FixedTimeStepping</type>
<t_initial>0</t_initial>
<t_end>1</t_end>
<!-- purely elastic problem: time scale irrelevant -->
<timesteps>
<pair>
<repeat>100</repeat>
<delta_t>0.01</delta_t>
<!-- i.e.: repeat a hundred times a time step with size of delta_t -->
</pair>
</timesteps>
</time_stepping>
</process>
</processes>
<output>
<type>VTK</type>
<prefix>OGSoutput_basin{:process_id}</prefix>
<timesteps>
<pair>
<repeat>100</repeat>
<each_steps>1</each_steps>
<!-- i.e.: do a hundred times the output at each step -->
</pair>
</timesteps>
<variables>
<variable>displacement</variable>
<variable>sigma</variable>
</variables>
<suffix>_ts_{:timestep}_t_{:time}</suffix>
</output>
</time_loop>
<parameters>
<!-- Material parameters -->
<parameter>
<name>YoungModulus</name>
<type>Group</type>
<group_id_property>MaterialIDs</group_id_property>
<index_values>
<index>0</index>
<value>60e6</value> <!--Pa-->
</index_values>
<index_values>
<index>1</index>
<value>60e7</value> <!--Pa-->
</index_values>
<index_values>
<index>2</index>
<value>60e6</value> <!--Pa-->
</index_values>
<index_values>
<index>3</index>
<value>60e9</value> <!--Pa-->
</index_values>
</parameter>
<parameter>
<name>PoissonRatio</name>
<type>Group</type>
<group_id_property>MaterialIDs</group_id_property>
<index_values>
<index>0</index>
<value>0.25</value>
</index_values>
<index_values>
<index>1</index>
<value>0.25</value>
</index_values>
<index_values>
<index>2</index>
<value>0.25</value>
</index_values>
<index_values>
<index>3</index>
<value>0.45</value>
</index_values>
</parameter>
<parameter>
<name>rho_sr</name>
<type>Constant</type>
<value>1</value>
</parameter>
<!-- Initial and boundary values -->
<parameter>
<name>displacement0</name>
<type>Constant</type>
<values>0 0</values>
</parameter>
<parameter>
<name>dirichlet0</name>
<type>Constant</type>
<value>0</value>
</parameter>
</parameters>
<process_variables>
<process_variable>
<name>displacement</name>
<components>2</components>
<order>1</order>
<initial_condition>displacement0</initial_condition>
<boundary_conditions>
<boundary_condition>
<mesh>mesh_basin_physical_group_Left</mesh>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
<boundary_condition>
<mesh>mesh_basin_physical_group_Right</mesh>
<type>Dirichlet</type>
<component>0</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
<boundary_condition>
<mesh>mesh_basin_physical_group_Top</mesh>
<type>Python</type>
<component>1</component>
<bc_object>bc_y</bc_object> <!--the python object is created using timeBCs_glacier.py (see L. 10)-->
<flush_stdout>true</flush_stdout> <!-- for debugging: false -->
</boundary_condition>
<boundary_condition>
<mesh>mesh_basin_physical_group_Bottom</mesh>
<type>Dirichlet</type>
<component>1</component>
<parameter>dirichlet0</parameter>
</boundary_condition>
</boundary_conditions>
</process_variable>
</process_variables>
<nonlinear_solvers>
<nonlinear_solver>
<name>basic_newton</name>
<type>Newton</type>
<max_iter>20</max_iter>
<linear_solver>general_linear_solver</linear_solver>
</nonlinear_solver>
</nonlinear_solvers>
<linear_solvers>
<linear_solver>
<name>general_linear_solver</name>
<eigen>
<solver_type>SparseLU</solver_type>
<scaling>true</scaling>
</eigen>
</linear_solver>
</linear_solvers>
</OpenGeoSysProject>
web/content/docs/tools/getting-started/video-tutorial/glacierclass.py 0 → 100644
+ 70
− 0
View file @ 99cd534d
# model of the evolving glacier extensions (length and height)
# parameterized, independent of concrete geometry
import numpy as np
import matplotlib.pyplot as plt
from math import pi, sin, cos, sinh, cosh, sqrt
gravity = 9.81 # m/s²
class glacier:
def __init__(self, L_dom, L_max, H_max, x_0, t_0, t_1):
self.rho_ice = 900 # kg/m³
self.L_dom = L_dom
self.L_max = L_max
self.H_max = H_max
self.x_0 = x_0
self.t_0 = t_0
self.t_1 = t_1
def normalstress(self, x, t):
return -self.rho_ice * gravity * self.local_height(x, t)
# analytical function for the glacier's shape
def local_height(self, x, t):
l = self.length(t)
if l == 0:
return 0 * x
else:
xi = (x - self.x_0) / l
xi = np.array(xi)
xi[xi > 1] = 1.0
return self.height(t) * np.sqrt(1 - xi ** 1)
def height(self, t):
return self.H_max * (t - self.t_0) / self.t_1
def length(self, t):
return self.L_max * (t - self.t_0) / self.t_1
def printMaxLoads(self):
print("Maximal normal stress due to glacier load: ")
print(self.normalstress(0, self.t_1) / 1e6, "MPa")
def plotEvolution(self):
tRange = np.linspace(self.t_0, self.t_1, 11)
xRange = np.linspace(self.x_0, self.x_0 + self.L_dom, 500)
yRangeBefore = 0
fig, ax = plt.subplots()
ax.set_title("Glacier evolution")
for t in tRange:
yRange = self.local_height(xRange, t)
ax.plot(xRange, yRange, label=t)
ax.fill_between(xRange, yRangeBefore, yRange)
yRangeBefore = yRange
ax.set_xlabel("$x$ / m")
ax.set_ylabel("height")
ax.grid()
fig.legend()
fig.savefig("glacier.pdf")
plt.show()
fig, ax = plt.subplots()
ax.plot(tRange, self.height(tRange))
ax.set_xlabel("$t$ / a")
ax.set_ylabel("height")
ax.grid()
web/content/docs/tools/getting-started/video-tutorial/history.sh 0 → 100644
+ 11
− 0
View file @ 99cd534d
#create the mesh via python control for gmsh
python3 mesh_basin.py
#transform the mesh and extract boundary meshes
python3 msh2vtu.py mesh_basin.msh --ogs --rdcd
#run OpenGeoSys (with the debug level information)
ogsr -l debug OGSinput_basin.prj
#do the postprocessing with ParaView
paraview OGSoutput_basin0.pvd
web/content/docs/tools/getting-started/video-tutorial/index.md 0 → 100644
+ 48
− 0
View file @ 99cd534d
+++
date = "2018-02-27T11:00:13+01:00"
title = "Video Tutorial"
author = "Dominik Kern"
weight = 105
[menu.docs]
name = "Tools & Workflows"
identifier = "tools"
weight = 4
post = "Helpful tools for pre- and postprocessing as well as complete model setup workflows."
[menu.tools]
parent = "getting-started"
+++
In this [three-part tutorial](https://www.youtube.com/watch?v=BULunRJQRJ0&list=PLU_clTnZqNAeOXENl79kQwn0pgHGittX1) we look at the mechanical deformation of a sedimentary basin under an advancing glacier. In detail we show you a possibility how to
* create a mesh with [Gmsh](http://gmsh.info/),
* implement space- and time-dependent boundary conditions,
* write an OGS input file,
* visualize the results with [ParaView](https://www.paraview.org/).
As a common programming language we use [Python](https://www.python.org).
## Part 1: Preprocessing
{{< youtube BULunRJQRJ0 >}}
## Part 2: Solving
{{< youtube GL5sugIyHEk >}}
## Part 3: Postprocessing
{{< youtube bkmubABAA_s >}}
## Supplementary material:
* [mesh_basin.msh](mesh_basin.msh)
* [mesh_basin.py](mesh_basin.py)
* [OGSinput_basin.prj](OGSinput_basin.prj)
* [timeBCs_glacier.py](timeBCs_glacier.py)
* [glacierclass.py](glacierclass.py)
* [history.sh](history.sh)
* [msh2vtu-project on GitHub](https://github.com/dominik-kern/msh2vtu)
* *Transient hydrodynamics within intercratonic sedimentary basins during glacial cycles*, V. F. Bense M. A. Person, DOI: [10.1029/2007JF000969](https://doi.org/10.1029/2007JF000969)
web/content/docs/tools/getting-started/video-tutorial/mesh_basin.msh 0 → 100644
+ 5648
− 0
View file @ 99cd534d
This diff is collapsed.
web/content/docs/tools/getting-started/video-tutorial/mesh_basin.py 0 → 100644
+ 226
− 0
View file @ 99cd534d
# ------------------------------------------------------------------------------
#
# Gmsh Python tutorial for meshing a layered sediment basin stratigraphy
#
# Geometry basics, elementary entities, physical groups
#
# ------------------------------------------------------------------------------
# The Python API is entirely defined in the `gmsh.py' module (which contains the
# full documentation of all the functions in the API):
import numpy
import gmsh
# Before using any functions in the Python API, Gmsh must be initialized:
gmsh.initialize()
# By default Gmsh will not print out any messages: in order to output messages
# on the terminal, just set the "General.Terminal" option to 1:
gmsh.option.setNumber("General.Terminal", 1)
# Next we add a new model named "t1" (if gmsh.model.add() is not called a new
# unnamed model will be created on the fly, if necessary):
gmsh.model.add("basinMesh")
# The Python API provides direct access to each supported geometry kernel. The
# built-in kernel is used in this first tutorial: the corresponding API
# functions have the `gmsh.model.geo' prefix.
# The first type of `elementary entity' in Gmsh is a `Point'. To create a point
# with the built-in geometry kernel, the Python API function is
# gmsh.model.geo.addPoint():
# - the first 3 arguments are the point coordinates (x, y, z)
# - the next (optional) argument is the target mesh size (the "characteristic
# length") close to the point
# - the last (optional) argument is the point tag (a stricly positive integer
# that uniquely identifies the point)
# The distribution of the mesh element sizes will be obtained by interpolation
# of these characteristic lengths throughout the geometry. Another method to
# specify characteristic lengths is to use general mesh size Fields (see
# `t10.py'). A particular case is the use of a background mesh (see `t7.py').
#
# If no target mesh size of provided, a default uniform coarse size will be used
# for the model, based on the overall model size.
# Units: meter
km = 1e3 # km in m
w = 120 * km # width
d = 7.5 * km # depth
lc = 1e3 # m
# We can then define some additional points. All points should have different tags:
gmsh.model.geo.addPoint(w / 2, 0, 0, lc, 1)
gmsh.model.geo.addPoint(-w / 2, 0, 0, lc, 2)
gmsh.model.geo.addPoint(-w / 2, -d, 0, lc, 3)
gmsh.model.geo.addPoint(w / 2, -d, 0, lc, 4)
# If the tag is not provided explicitly, a new tag is automatically created, and
# returned by the function:
pLay1centr = gmsh.model.geo.addPoint(0, -d / 6, 0, lc)
pLay2centr = gmsh.model.geo.addPoint(0, -d / 3, 0, lc)
pLay3centr = gmsh.model.geo.addPoint(0, -d / 2, 0, lc)
# Auxiliary points for the sediment layers
pLay1north = gmsh.model.geo.addPoint(-w / 8, 0, 0, lc)
pLay1south = gmsh.model.geo.addPoint(w / 8, 0, 0, lc)
pLay2north = gmsh.model.geo.addPoint(-w / 6, 0, 0, lc)
pLay2south = gmsh.model.geo.addPoint(w / 6, 0, 0, lc)
pLay3north = gmsh.model.geo.addPoint(-w / 4, 0, 0, lc)
pLay3south = gmsh.model.geo.addPoint(w / 4, 0, 0, lc)
# Curves are Gmsh's second type of elementery entities, and, amongst curves,
# straight lines are the simplest. The API to create straight line segments with
# the built-in kernel follows the same conventions: the first 2 arguments are
# point tags (the start and end points of the line), and the last (optional one)
# is the line tag.
#
# Note that curve tags are separate from point tags - hence we can reuse tag `1'
# for our first curve. And as a general rule, elementary entity tags in Gmsh
# have to be unique per geometrical dimension.
# gmsh.model.geo.addLine(1, 2, 1) # (start, end, tag)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)
gmsh.model.geo.addLine(pLay3south, pLay2south, 5)
gmsh.model.geo.addLine(pLay2south, pLay1south, 6)
gmsh.model.geo.addLine(pLay1south, pLay1north, 7)
gmsh.model.geo.addLine(pLay1north, pLay2north, 8)
gmsh.model.geo.addLine(pLay2north, pLay3north, 9)
gmsh.model.geo.addLine(1, pLay3south, 10)
gmsh.model.geo.addLine(pLay3north, 2, 14)
gmsh.model.geo.addSpline([pLay1south, pLay1centr, pLay1north], 11)
gmsh.model.geo.addSpline([pLay2south, pLay2centr, pLay2north], 12)
gmsh.model.geo.addSpline([pLay3south, pLay3centr, pLay3north], 13)
# The third elementary entity is the surface. In order to define a surface
# from the curves defined above, a curve loop has first to be defined.
# A curve loop is defined by an ordered list of connected curves,
# a sign being associated with each curve (depending on the orientation of the
# curve to form a loop). The API function to create curve loops takes a list
# of integers as first argument, and the curve loop tag (which must be unique
# amongst curve loops) as the second (optional) argument:
# e.g. gmsh.model.geo.addCurveLoop([4, 1, 2, 3], 1)
gmsh.model.geo.addCurveLoop([7, -11], 1)
gmsh.model.geo.addCurveLoop([6, 11, 8, -12], 2)
gmsh.model.geo.addCurveLoop([5, 12, 9, -13], 3)
gmsh.model.geo.addCurveLoop([2, 3, 4, 10, 13, 14], 4)
# Add plane surfaces defined by one or more curve loops. The first curve
# loop defines the exterior contour; additional curve loop define holes.
# (only one here, representing the external contour, since there are no holes
# --see `t4.py' for an example of a surface with a hole):
for l in range(1, 5):
gmsh.model.geo.addPlaneSurface([l], l)
# At this level, Gmsh knows everything to display the surfaces and
# to mesh it. An optional step is needed if we want to group elementary
# geometrical entities into more meaningful groups, e.g. to define some
# mathematical ("domain", "boundary"), functional ("left wing", "fuselage") or
# material ("steel", "carbon") properties.
#
# Such groups are called "Physical Groups" in Gmsh. By default, if physical
# groups are defined, Gmsh will export in output files only mesh elements that
# belong to at least one physical group. (To force Gmsh to save all elements,
# whether they belong to physical groups or not, set the `Mesh.SaveAll' option
# to 1.) Physical groups are also identified by tags, i.e. stricly positive
# integers, that should be unique per dimension (0D, 1D, 2D or 3D). Physical
# groups can also be given names.
#
# Here we define physical curves that groups the left, bottom, top and right
# curves in single groups for the one-dimensional boundary meshes
dim = 1
North = gmsh.model.addPhysicalGroup(dim, [2])
gmsh.model.setPhysicalName(dim, North, "Left")
Base = gmsh.model.addPhysicalGroup(dim, [3])
gmsh.model.setPhysicalName(dim, Base, "Bottom")
South = gmsh.model.addPhysicalGroup(dim, [4])
gmsh.model.setPhysicalName(dim, South, "Right")
Ground = gmsh.model.addPhysicalGroup(dim, [10, 5, 6, 7, 8, 9, 14])
gmsh.model.setPhysicalName(dim, Ground, "Top")
# ... and we define physical surfaces for the different
# material domains containing the corresponding geometrical surfaces:
dim = 2
Lay1 = gmsh.model.addPhysicalGroup(dim, [1])
gmsh.model.setPhysicalName(dim, Lay1, "SedimentLayer1")
Lay2 = gmsh.model.addPhysicalGroup(dim, [2])
gmsh.model.setPhysicalName(dim, Lay2, "SedimentLayer2")
Lay3 = gmsh.model.addPhysicalGroup(dim, [3])
gmsh.model.setPhysicalName(dim, Lay3, "SedimentLayer3")
Bed = gmsh.model.addPhysicalGroup(dim, [4])
gmsh.model.setPhysicalName(dim, Bed, "RockBed")
# Before it can be meshed, the internal CAD representation must be synchronized
# with the Gmsh model, which will create the relevant Gmsh data structures. This
# is achieved by the gmsh.model.geo.synchronize() API call for the built-in
# geometry kernel. Synchronizations can be called at any time, but they involve
# a non trivial amount of processing; so while you could synchronize the
# internal CAD data after every CAD command, it is usually better to minimize
# the number of synchronization points.
gmsh.model.geo.synchronize()
# We can then generate a 2D mesh...
gmsh.model.mesh.generate(2)
# ... and save it to disk: msh=current format, msh22 older format
gmsh.write("mesh_basin.msh")
# Remember that by default, if physical groups are defined, Gmsh will export in
# the output mesh file only those elements that belong to at least one physical
# group. To force Gmsh to save all elements, you can use
#
# gmsh.option.setNumber("Mesh.SaveAll", 1)
# By default, Gmsh saves meshes in the latest version of the Gmsh mesh file
# format (the `MSH' format). You can save meshes in other mesh formats by
# specifying a filename with a different extension. For example
# gmsh.write("basinMesh.pdf") only accessible via GUI
# gmsh.write("basinMesh.pvtu")
# gmsh.write("basinMesh.vtk")
# gmsh.write("basinMesh.unv")
# gmsh.write("basinMesh.stl")
# will save the mesh in the vtk, unv and stl format.
# To visualize the model we can run the graphical user interface with:
gmsh.fltk.run()
# Note that starting with Gmsh 3.0, models can be built using other geometry
# kernels than the default "built-in" kernel. To use the OpenCASCADE geometry
# kernel instead of the built-in kernel, you should use the functions with the
# `gmsh.model.occ' prefix.
#
# Different geometry kernels have different features. With OpenCASCADE, instead
# of defining the surface by successively defining 4 points, 4 curves and 1
# curve loop, one can define the rectangular surface directly with
#
# gmsh.model.occ.addRectangle(.2, 0, 0, .1, .3)
#
# Boolean operation using OpenCascade module (occ, see below)
# gmsh.model.occ.cut([(3, 1)], [(3, 2)], 3)
#
# After synchronization with the Gmsh model with
#
# gmsh.model.occ.synchronize()
#
# the underlying curves and points could be accessed with
# gmsh.model.getBoundary().
#
# See e.g. `t16.py', `t18.py', `t19.py' or `t20.py' for complete examples based
# on OpenCASCADE, and `demos/api' for more.
# This should be called when you are done using the Gmsh Python API:
gmsh.finalize()
web/content/docs/tools/getting-started/video-tutorial/timeBCs_glacier.py 0 → 100644
+ 35
− 0
View file @ 99cd534d
import OpenGeoSys
import glacierclass as glc
L_dom = 120000 # m
L_max = 0.7 * L_dom # m
H_max = 200 # m
x_0 = -0.5 * L_dom # m
t_0 = 0.00 # s
t_1 = 1.0000 # s
class BC_Y(OpenGeoSys.BoundaryCondition):
def __init__(self, L_dom, L_max, H_max, x_0, t_0, t_1):
super(BC_Y, self).__init__()
# instantiate the glacier member object
self.glacier = glc.glacier(L_dom, L_max, H_max, x_0, t_0, t_1)
self.glacier.printMaxLoads()
self.glacier.plotEvolution()
def getFlux(
self, t, coords, primary_vars
): # here Neumann BC: flux of linear momentum
x, y, z = coords
# print("x = ", x)
# print("t = ", t)
# print("s = ", self.glacier.normalstress_glacier(x,t))
value = self.glacier.normalstress(x, t)
derivative = [0.0, 0.0]
return (True, value, derivative)
# instantiate the BC objects used by OpenGeoSys
bc_y = BC_Y(L_dom, L_max, H_max, x_0, t_0, t_1)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment