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
            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 ...
            where you can change the sub-area without stepping time:
INSPECT>i #i1,#i2,#j1,#j2 [c|n]
            To print the whole area:
            (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 ...
            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   - 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.

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

bold: fixed text, italic: user parameters, # indicates that a numerical value is expected, [ ]: optional.

 -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).
 -Od1,d2      - output devices for u-array,v-array; see above,
                "Output to other xterm..."
output sub-array (default: entire array)
[#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 (!))
[#a]     - output type "Stress", Mohr-Coulomb stress margin instead of σxx,
Specify amplitude factor a (default: 1e-9),
witches Mohr-Coulomb option on.
                A non-numeric value will switch off this option (i.e. will
σxx again).
 -oc[#a]      - Character translation of displacement or stress with amplification a.
  -on[#a]      - Numeric output of displacement with amplification a.

 -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.
(no parameters) switches off.
 -mcON          Switch on with default/current parameters.

 -li#r,#t,#s  - lithostatic conditions: σyy = (s/dxg·ρ·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.

Interior stresses on a cube around i+1/2, j+1/2

              - Normal surface stress =
σyy (upper boundary) raised cosine
                amplitude f from location i1 to i2, duration 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


Start two character display windows (xterm):
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 boundar
y 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
-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      
  9                                           AAAA       aaaa       
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