2dwavet - a test program for Finite Difference
code
2-dimensional seismic wave equation (and with friction:
elastostatic solver *))
*) incomplete solutions since the wave equation misses nontrivial
solutions of the stress-free case. Keyword: Airy stress function.
Without command line parameters, the solver processes
model size:
11 x 11
vp:
8 km/s
Poisson: 0.25
sub-crit: 0.8
seed: v=1 at
5,5 at first time step
boundary c.: U: free, L: symm, R and B: reflective
damping: No
Wave propagation problems can be solved (for that the
subroutines sw2d and bcw2d should be developed further) as well as
elasto-static cases. For the latter, the damping coefficient should
be set to a value -0.5 or so in order to allow quick propagation of
the perturbation, avoid reflections, and reach convergence in
reasonable time.
A "viscous" problem can be solved for an elastostatic problem
with the relaxation method. This variant is selected with the -v
option.
Output and interaction:
Normally, U and V arrays are written to STDOUT at every step.
At the prompter ...
NEXT
#>
A few commands at
position 1 in the response line: i f
See INSPECT>.
To
stop
: Enter any other non-numeric character + <CR>
To continue : Press <CR>
To skip n
prints : Enter a number n + <CR>
To bypass
prompting n times: NEXT #> +A#n
To wake up graphics: NEXT #> +G [gr-options]
The +-specifications are
options; they need at least one blank
at start; they can be appended to all responses at NEXT
#>
See Graphics control
below.
To select a sub-area for inspection:
NEXT #>i
#i1,#i2,#j1,#j2 [c|n]
where c requests character form and n
numeric
You will get a prompter ...
INSPECT>
where you can change the sub-area without stepping time:
INSPECT>i #i1,#i2,#j1,#j2 [c|n]
To print the whole area:
INSPECT>f
(works also at the NEXT> prompter).
To continue with time-stepping, press <CR>
NEXT #>
To enter an option (one of the command line options)
code it a the NEXT #> prompter
You will receive the prompter ...
OPTION>
where you can enter additional options of the command-line category.
To continue with
time-stepping, press
<CR>
Always open a second xterm (diagnostic output)!
Dialog is kept on STDIN
STDOUT using
ncurses.
Graphics
control at the NEXT #> prompter:
The following options can be entered in combination with
other entries if they are appended to the end. At least one leading
blank is needed if no other option is entered. Options must be
separated with blanks.
By default, graphics is shown once after init/revive and must be
revived at any time desired.
It can be kept live for a period of time.
The other options are remembered.
Options can be combined up to a length of 63 characters.
Palette file default =
/home/hgs/Fdiff/palette
Palette files must be produced with the GMT makecpt command
You may consider
ln -s file.cpt /home/hgs/Fdiff/palette
+G[#n]
- revive graphics and keep alive n times.
-G
- skip graphics
+GR[file] - palette
refresh (with new file)
+U
-U - show U-displacement (don't show)
+V
-V - show V-displacement
+SXX -SXX
- show stress_xx
+SYY -SYY
- show stress_yy
+SXY -SXY
- show stress_xy
+M -M
- Switch Mohr-Coulomb analysis on/off. Off: the tectonic stress is
shown if +SXX
+TS -TS
- Switch Tectonic stress option on/off. Off: The deformation
stress is shown if +SXX
+SXX +MC - (together)
Show Mohr-Coulomb stress margin instead of stress_xx
(Why together? +MC needs the stress field, and we save MC-stress
in the SXX array)
(You would certainly request +SXX +SYY +MC)
+SXX -M +TS -
msc=#v - scale
Mohr-Coulomb stress margin with v
(no blanks before and after `=´ here and in the following
instructions)
ssc=#v - scale stress
with v
dsc=#v - scale
displacement with v
Output
to other xterm devices:
Start a third xterm. Determine devices with ps.
Assume you obtain pts/2 and pts/3
Open one device (unit number used: 11):
2dwavet -O/dev/pts/2
Open two devices (unit numbers used: 11, 12):
2dwavet
-O/dev/pts/2,/dev/pts/3
Short form (numeric arguments): -O2,-1 (one device)
-O2,3 (two)
The short form supposes /dev/pts/
as the device directory.
~/bin/2dww and
~/loxterms:
Collects /dev information about recently opened xterms
Starting a bash shell with $NOTEDEV = YES in the environment will
write a line to ~/loxterm,
e.g. /dev/pts/0
The shell script 2dww
uses this feature to open two xterms and collect the
device names in /home/hgs/loxterms
If the file exists, 2dwavet
will try to use the information to prepare
the two xterms for output. If that succeeds, you would avoid the -O
option.
If 2dwavet
fails, issue
rm ~/loxterms; 2dww
OPTIONS:
bold: fixed text, italic: user parameters, # indicates that a
numerical value is expected, [ ]:
optional.
Model:
-A#m,#n
- model size (default 11,11)
-N#n
- Number of time steps
-s#s
- Subcriticality (default: 0.8)
-DX#dx
- Grid constant, used for lithostatic stress and for
model dimensionality (default: 1.0)
-PG#g
- Gravity acceleration (litho stress)
-PS#s
- Stress scale. Default 1e9, i.e. stresses are shown
in units of Gpa. Where stresses are dealt with in
the diff equations,
physical value of displacement difference =
numeric value of stress * s
* dx /(λ+2μ)
(options -F -Fs and subroutines fw2d_ssyy
and fw2d_stensor).
Output:
-Od1,d2
- output devices for u-array,v-array; see above,
"Output to other xterm..."
-OR#i1,#i2,#j1,#j2
- output sub-array (default: entire array)
-OD[#a] -
output type "Displacements" U to device d1, V to d2
with amplification a
(default: 1.0)
-OS[#a]
- output type "Stress", ix·σxx + iy·σyy to device d1 and
σxy to d2. Specify amplitude factor a (default: 1e-9)
Is neutral to the Mohr-Coulomb and Tectonic stress options.
ix and iy are set with the -OSC option.
-OSC#ix,#iy
- select the combination of diagonal elements for -OS
-OS... implies always
character translation.
ix and iy are integers, preferably 0
and 1 (default: 1,1 (!))
-OMC[#a]
- output type "Stress", Mohr-Coulomb stress margin instead of σxx,
however. Specify
amplitude factor a (default:
1e-9),
switches Mohr-Coulomb
option on.
A non-numeric value will switch off this option (i.e. will
output σxx again).
-oc[#a] - Character
translation of displacement or stress with amplification a.
-on[#a] - Numeric output
of displacement with amplification a.
Material:
-vp#v
- compressional wave velocity [m/s]
-p#p | -mu#m -
Poisson constant (default: 0.25) | Shear modulus [Pa] (!)
-rho#r
- density ρ (default: 3000 kg/m3).
-d#d
- wave problem, damping viscosity d, a negative value !
-v#d
- static problem solved with relaxation.
Damping viscosity d,
must be greater than +1 !
-mc#f[,#t] - Mohr-Coulomb law
with friction f and
tension t.
Switches Mohr-Coulomb option on.
-mc
(no parameters) switches off.
-mcON
Switch on with default/current parameters.
-li#r,#t,#s -
lithostatic conditions: σyy
= (s/dx)·g·ρ·y,
σxx = r·σyy
+ (s/dx)·t
where y = (j-1/2)·dx.
Defaults: t = 0 [Gpa]
(!), s = 1000.
Switches Tectonic option on.
-li
(no parameters) switches Tectonic option off.
-liON
Switch on with default/current parameters.
Boundary conditions:
The top boundary is always stress-free
The right and bottom boundaries are reflective by default
(aux column/row = 0)
The left boundary causes symmetry by default.
-BR[b|l|r]
- radiating boundary bottom, left, right
-BF[u|b|r]
- xy and yy free boundary up, bottom, right
-B0l
- left boundary reflective, not symmetric
-BBB[#d] - Buoyant
bottom boundary, density jump d
[kg/m3], default 500.
Starting conditions:
-S
- use the default starting seed u(5,5)=1.0
-Su#i,#j,#u - a starting ("seed") value for the
u-array
-Sv#i,#j,#v
- a starting ("seed") value for the v-array
-S[u|v] can be repeated 16 times.
-Fi#i,#j,#sxx,#syy,#sxy
- Interior stresses on a cube around i+1/2, j+1/2
-Fn#i1,i2,#syy[,#n]
- Normal surface stress = σyy (upper boundary) raised cosine
amplitude f from
location i1 to i2, duration n.
-Fs#i1,i2[,#sxy[,#n]]
- Shear force couple upper and lower boundary.
Default value is computed from the normal force so that
a torque balance is achieved.
-g#v
- Balance vertical stress with body force, scale with factor v
(default 1.0)
-Fd#n - Force
duration
EXAMPLES:
Start two character display windows (xterm):
2dww
An elasto-static problem with buoyant bottom
boundary:
2dwavet -oc100. -Fn40,60,1.,1000 -A100,10 -DX1000. -d-0.3
-BBB500. -BRr
(Remember: the solution to the wave equation allows
static displacements,
solutions to the homogeneous Airy differential equations, so the
solution
of this program is incomplete)
This on is with only free boundaries. It does not
balance torque
2dwavet -oc1000. -Fn40,60,1.,1000 -A100,10 -DX1000. -d-0.3
-s0.3 -BFrbu
This one balances (tries to balance) torque:
2dwavet -oc1000. -Fn40,60,1.,1000 -Fs3,5 -A100,10 -DX1000.
-d-0.3 -s0.3 -BFrbu
(note that we don't specify the traction magnitude
in -Fs )
Result (at 81 converged to the end; note the
vertical shift at the left boundary and the warp to the right next
to it - could be a bug):
U: 41 2
1 @
6666@6662cjicDEBaa
ACJOk@@@@@@@@@@
6666666666Kojca
2 @
66666666VelhaFEAaa
ACFP66666666666
@@@@@@@@@@@pfca
3 @ 66666666PBbc
CCA
ABFO66666666666
@@@@@@@@@@@ofba
4 6 @@66666YJFDab
AA
ABEKX6666666666
@@@@@@@@@@xkeba
5 6 @@o666TEaaBB
a
AACFMW66666666Q
q@@@@@@@@wmfcaa
6 @ 66O@@@teAAbb
A
AAACEINSVXXVRMH
gmrvxxvsniecaaa
7 @ 66@@@@@yjfdAB
aa
AABDFHIIIGHEC
ceghiiigfdbaa
8 6 @@@@@@@@pbBC
cca
AABBCDEEFFEDCB
bcdeffeedcbbaa
9 6
@@@@@@@@vELGAfeaAA
AABBBCCCCCCCBBA
abbcccccccbbbaa
10 6
@@@@6@@@=CJICdebAA
AAABBCCDEEEEEDCA
acdeeeeedccbbaaa
>>>
.........10........20........30........40........50........60........70........80........90.......100
U: 81 3
1 @
6666@666ZIDAA
AADJPk@@@@@@@@@@
6666666666Kpjdaa
2 @
66666666VICAA
ACHP66666666666
@@@@@@@@@@@pgca
3 @
66666666RHCA
ACFP66666666666
@@@@@@@@@@@pfca
4 6
@@666665LEBA
ABELY6666666666
@@@@@@@@@@yleba
5 6
@@v666YKDBA
ABCHNY66666666Q
q@@@@@@@@yngcba
6 @
66V@@@ykdba
AABDGLRW133ZVPG
hpvz&&+wrlhdbaa
7 @
66@@@@@#leba
AABDFGKLMMLJHD
dgjlmmlkhfdbaa
8 6
@@@@@@@@rgca
AABCDEFFFEDCB
bcdefffedcbaa
9 6
@@@@@@@@vicaa
AABBCCDDCCBA
abccddccbbaa
10 6
@@@@6@@@zidaa
AABCDEFHHFEDB
bdefggfedcbaa
>>>
.........10........20........30........40........50........60........70........80........90.......100
V: 41 2
1 666@@@@666MFbge
BB
aadmm666666666666666666666mmdaa
2 666@@@@66L
cebAA
bga666666666666666666666agb
3
6666@@@@hgjd
aa
AK666666666666666666666KA
4
6666@@@@@ukbBAaba
ADM666666666666666666666MDA
5 6666@@@@@=jDF
cba
ABDKY6666666666666666666YKDBA
6 6666@@@@@=jDF
cba
ACHO1666666666666666661OHCA
7
6666@@@@@ukbBAaba
ABDHMT366666666666663TMHDBA
8
6666@@@@hgjd
aa
AABCEGKOSWZ34543ZWSOKGECBAA
9 666@@@@66L
cebAA
AAABCEFGJLMNNNMLJGFECBAAA
10 666@@@@666MFbge
BB
AAAAAAABCEFGJLMNNNMLJGFECBAAAAAAA
>>>
.........10........20........30........40........50........60........70........80........90.......100
V: 81 3
1
666@@@@664IDBA
abemm666666666666666666666mmeba
2
666@@@@66IA
bga666666666666666666666agb
3
6666@@@@uhdba
AL666666666666666666666LA
4
6666@@@@@thcaa
ADM666666666666666666666MDA
5
6666@@@@@yjdba
ABDJY6666666666666666666YJDBA
6
6666@@@@@yjdba
AACHNZ66666666666666666ZNHCAA
7
6666@@@@@thcaa
ABDGNU566666666666665UNGDBA
8
6666@@@@uhdba
AABDHKPUZ566666665ZUPKHDBAA
9
666@@@@66IA
ABCEHJMPSVXYYYXVSPMJHECBA
10
666@@@@664IDBA
AABCEHJMPSUWYYYWUSPMJHECBAA
>>>
.........10........20........30........40........50........60........70........80........90.......100
We could try with a normal traction that varies linearly along
the left boundary from tension (y=1) to compression (y=10). We
cannot do it with symmetry switched on and surface shear.
Then, a normal bottom traction near x=2 would do to balance
the two surface normal loads. Could be interesting.
New call: undamped wave problem
2dwavet -oc100000. -Fn40,60,1.,2 -A100,50 -BFu -BRr
Old call: A wave problem
2dwavet -O2,3 -oc50.
-F40,60,1.,100 -d-0.1 -A100,50
U:
13 1
1
abbcccbbaa AABBCCCBBA
2
AABCCCCBBA abbccccbaa
3
ABBBBBBAA aabbbbbba
4
aaaaaaaaa
AAAAAAAAA
5
aabbbbcbbaa AABBCBBBBAA
6
aabbbbbbaa AABBBBBBAA
7
aaaaaaa AAAAAAA
8
9
AAAA
aaaa
10
11
V:
13 1
1
BHOZ6666666666666ZOHB
2
ACHLQW366666663WQLHCA
3
ACEHJLNOPQPONLJHECA
4
AABBCCDDDDDCCBBAA
5
aaaabbbbbbbbbbbaaaa
6
aaabbcddeefffeeddcbbaaa
7
abcdeghijkkkjihgedcba
8
aabcdeeefeeedcbaa
9
AAAABBBBBBBBBBBAAAA
10
AAABBBCCCCCBBBAAA
11
AAAAAAAAA
12
13
.bye