Some Notes about Grid Generation
OVERVIEW
Philippe F. has written code to generate a grid. This code runs as part of the
interactive pre-processing part of the thermo-mechanical model.
The grid generation code allows the user to define the different sections of the
grid by defining the characteristics of the "interface" between each part of the
grid. The code then uses these definitions to interpolate and create a grid.
Figure 1. Example of a desired grid.
(have one or two examples, covering all the types of interface)
INPUT
The grid generation code asks for the following parameters:
hx width of plotting window
hz height of plotting window
hx, hz are the width and height of plotting window (in the same units
as used for the grid - normally metres)
(0,0) is the assumed origin at upper left
NOTE: hx should be the same width that you want your grid to be.
So if you want a grid 1000 km wide, use hx of 1000000.
igy number of nodes in the horizontal direction, for the thermal grid
igz number of nodes in the vertical direction, for the thermal grid
You are then asked for the number of lines. This is the number of
interfaces in your grid. In the example below, there are 4 of these
interfaces.
For each interface, you need to enter:
* The row number in the grid for this interface
* What type of interface:
a) line-point
b) translate previous interface
c) create a constant thickness slab
d) line-circle-line interface
e) circular transition from a horizontal to a horizontal
(e.g subduction) geometry
f) read-in interface
g) line segments
(these types of interfaces are explained below)
* Parameters needed to create the chosen type of interface (explained
below)
TYPES OF INTERFACES
a) line-point
This type of interface is a straight line. You define it by giving the
(x,z) coordinates of the first point on the line, followed by the
(x,z) coordinates of the second point on the line, followed by the number
of points on the line.
The number of points used determines how well the resulting grid lines
follow the desired interface. More points, better fit... but you must
keep in mind that there is a limit to the total number of points in
all the interfaces. (Looks like 1000 points at the moment - BL)
In general, you should put the endpoints of the line slightly outside the
desired grid area, so the code can extrapolate accurately at the edges of the
grid.
Example:
enter reference line npoints and xa za xb zb
0 0 605000 0 201
This defines a line from (0,0) to (0,605000) with 201 points on the line.
NOTE for this example: The prompt is wrong! The order of the input should
be: xa za xb zb npoints
where xa is the x coordinate of the first point of the segment
za is the z coordinate of the first point of the segment
xb is the x coordinate of the end point of the segment
zb is the z coordinate of the end point of the segment
npoints is the number of points on the line
The "first" point is the one with the smaller x value.
GENERAL NOTE: Philippe says, be careful that the total number of
points on all interface segments is not more than 1000.
b) translate previous interface
This creates an interface which is the same shape as the previous one, but
translated in the z direction.
[ Add more info here]
c) create a constant thickness slab
d) line-circle-line interface
(diagram)
Use this interface type when one or both of the lines are not horizontal. If
both lines are horizontal, use "circular transition from a horizontal to a
horizontal geometry".
e) circular transition from a horizontal to a horizontal (e.g. subduction)
geometry
(diagram)
NOTE:
For horizontal-circle-horizontal transition, can only specify horizontal
position of singularity, exit angle, and depth of slab. Can't specify the
delta x (horizontal distance spanned by the arc); have to let the program
place the top of the circle.
Example:
enter horizontal singularity position
300000
enter zleft zright
55588 20000 [use positive numbers; "z" here is not right; we're looking at
depth]
enter slab angle
20
enter slab length
500000 [use a big number to make sure it goes past edge of grid]
radius is 590111.1932229433 55588.0 20000.0
confirm :delta_z= 35588.0
top circle position 501829.7513867472 -20000.0
sla1 point 300000.0 -55588.0
enter angle left
0
500000.0 300000.0 3.14159 353553.625139613 -200000.0
-200000.0 -55588.0 300000.0 -55588.0
enter points on exit-slab ,circle,in-slab
101 31 101
enter position of top of in-slab
605000 -20000
f) read-in interface
[check the code about this... or ask Philippe]
==========================================================================
INPUT file notes:
in microfem_tuning_i, set ifpa = 1 if we want to generate result_track
(ifpa = i first pass)
[cobra]/disk4/bonny/ptt/c13q_v1_e1-1.75my_rb1 % ../CODE960916/xnew
1=1/3 formul 2 = 1/2 formul
1
1e+38= 1.000000000000000E+38
uncoupled run
TIME INITIALIZATION : .0 473041372198.4124
igridn 0
201 28 5628 301 28 8428 301 55
generate (1=interactive 2=file) or read (0) mesh
1
make thermal grid .choose here your hx hz
exemple 400000. 100000.
[hx, hz are width and height of plotting window (in the same units
as used for the grid - normally metres)
(0,0) is the assumed origin at upper left]
600000 300000
your z is exagerated (or equivalently your x
is shrunk by a factor 4.0
enter igy igz
[igy is number of nodes in horizontal direction in the thermal grid.
igz is number of nodes in vertical direction in the thermal grid. ]
201 55
io_open_xwindow: Opening X11 window
[ displays window with title Plot "E mesh 2" with a grid]
enter number of lines
[ number of lines = number of depth interfaces ]
[ "line" (ligne) = "row" ]
4
line number 1
enter line number in grid
1
new menu !!!
enter 0 for line-point
enter 1 to translate previous interface
enter 2 to create a constant thickness slab
enter 3 for line-circle-line interface
enter 5 for circular transition from
a horizontal to a horizontal (e.g subduction)
geometry
enter 4 for read-in interface
0
enter reference line npoints and xa za xb zb
[ Actually this is not the right order for input.
Enter xa za first (coordinates of first point of segment)
in order of increasing x
Enter xb zb next (coordinates of end point of segment)
Enter npoints (number of points on the segment)
NOTE: be careful that the total number of points on all interface
segments is not more than 1000.
]
0 0 605000 0 201
sub menu along interface ! : [in plot window, see interface drawn]
1= get points on screen
0= translate previous
-1 = add no points
[ these options give the ability to add further points to this interface
by hand]
-1
[ ... coordinates of points on interface ...]
line number 2
enter line number in grid
28
new menu !!!
enter 0 for line-point
enter 1 to translate previous interface
enter 2 to create a constant thickness slab
enter 3 for line-circle-line interface
enter 5 for circular transition from
a horizontal to a horizontal (e.g subduction)
geometry
enter 4 for read-in interface
5
enter horizontal singularity position
300000
enter zleft zright
55588 20000
enter slab angle
20
enter slab length
500000 [big number]
radius is 590111.1932229433 55588.0 20000.0
confirm :delta_z= 35588.0
top circle position 501829.7513867472 -20000.0
sla1 point 300000.0 -55588.0
enter angle left
0
500000.0 300000.0 3.14159 353553.625139613 -200000.0
-200000.0 -55588.0 300000.0 -55588.0
enter points on exit-slab ,circle,in-slab
101 31 101
enter position of top of in-slab
605000 -20000
[lots of coordinates...]
line number 3
enter line number in grid
43
new menu !!!
enter 0 for line-point
enter 1 to translate previous interface
enter 2 to create a constant thickness slab
enter 3 for line-circle-line interface
enter 5 for circular transition from
a horizontal to a horizontal (e.g subduction)
geometry
enter 4 for read-in interface
3
enter top circle position. radius, angle
500000 -21000 590111.19 20
enter slab length
500000
500000.0 298170.2497155633 3.14159 353553.625139613 -171676.111098546
-171676.111098546 -227597.932937419 298170.2497155633
-56587.9998056331
enter points on exit-slab ,circle,in-slab
50 30 50
enter position of top of in-slab
605000 -21000
[...lots of coordinates ...]
line number 4
enter line number in grid
55
new menu !!!
enter 0 for line-point
enter 1 to translate previous interface
enter 2 to create a constant thickness slab
enter 3 for line-circle-line interface
enter 5 for circular transition from
a horizontal to a horizontal (e.g subduction)
geometry
enter 4 for read-in interface
2
enter slab thickness
60000
[ points written out ]
z -interpolate
now plot
define line number that you consider as
the top of the descending slab
43
done
E2 grid
field range .0 .0
No variations in field
velocity range is : .0
No variation in velocity field
[ Plot window (xplot) showing thermal grid
Click on "change buttons"
Click on "save grid"
]
saving grid into file save_grid_i
enter row number of upper slab if relevant
43
done ! you can now rename the file
redefine grid E 1
E1 grid
field range .0 .0
No variations in field
velocity range is : .0
No variation in velocity field
[ displays E1 grid in xplot window ]
[ From here on, same as for SUBDUCT960826 instructions ]