Extended fixes to branch cdm_06262019

models/cdm/cdm/libcdm.py

@@ -219,7 +219,9 @@ def CDM(X, Y, X0, Y0, depth, ax, ay, az, omegaX, omegaY, omegaZ, opening, nu):
                            R3[2], R4[2]])
 
     if numpy.any(VertVec>0):
-        raise ValueError('Half-space solution: The CDM must be under the free surface!' +
+        # raise ValueError('Half-space solution: The CDM must be under the free surface!' +
+        #                  ' VertVec={}'.format(VertVec))
+        print('Half-space solution: The CDM must be under the free surface!' +
                          ' VertVec={}'.format(VertVec))
 
     if ax == 0 and ay == 0 and az == 0:

models/cdm/examples/cdm.pfg

@@ -1,12 +1,12 @@
 ;
 ; michael a.g. aïvázis
 ; orthologue
-; (c) 1998-2019 all rights reserved
+; (c) 1998-2018 all rights reserved
 ;
 
 cdm:
-    ; for running with MPI
-    ; shell = mpi.mpirun
+    ; shell
+    ; shell = mpi.shells.mpirun ; for running with mpi
 
     controller:
         ; saving the results of the calculation
@@ -17,7 +17,6 @@ cdm:
 
     ; job layout
     job:
-        hosts = 1       ; number of MPI hosts
         tasks = 1       ; number of tasks per host
         gpus = 0        ; number of gpus per task
         chains = 2**12  ; number of chains per task
@@ -46,54 +45,53 @@ cdm:
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (-5000, 5000)
+                    support = (-1000, 1000)
                 prior:
-                    support = (-5000, 5000)
+                    support = (-1000, 1000)
             depth:
                 count = 1
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (2000, 4000)
+                    support = (4000, 10000)
                 prior:
-                    support = (2000, 4000)
+                    support = (4000, 10000)
             opening:
                 count = 1
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (0, 1)
+                    support = (10, 80)
                 prior:
-                    support = (0, 1)
+                    support = (10, 80)
             a:
                 count = 3
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (0, 1)
+                    support = (100, 3000)
                 prior:
-                    support = (0, 1)
+                    support = (100, 3000)
             omega:
                 count = 3
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (0, 45)
+                    support = (-5, 5)
                 prior:
-                    support = (0, 45)
+                    support = (-5, 5)
             offsets:
-                count = 2
+                count = 1 
                 prep = uniform
                 prior = uniform
                 prep:
-                    support = (-0.1, 0.1)
+                    support = (-1.0, 1.0)
                 prior:
-                    support = (-0.1, 0.1)
+                    support = (-1.0, 1.0)
 
 
 ; for parallel runs
-mpi.shells.mpirun # cdm.shell:
-    hostfile = localhost
+mpi.shells.mpirun # altar.plexus.shell:
     extra = -mca btl self,tcp
 
 ; end of file

models/cdm/examples/synthetic/cdm.py

@@ -31,28 +31,28 @@ class CDM(altar.application, family="altar.applications.cdm"):
     y = altar.properties.float(default=0)
     y.doc = "the y coördinate of the CDM source"
 
-    d = altar.properties.float(default=3000)
+    d = altar.properties.float(default=7000)
     d.doc = "the depth of the CDM source"
 
-    aX = altar.properties.float(default=400)
+    aX = altar.properties.float(default=1000)
     aX.doc = "the x semi-axis length"
 
-    aY = altar.properties.float(default=450)
+    aY = altar.properties.float(default=1000)
     aY.doc = "the y semi-axis length"
 
-    aZ = altar.properties.float(default=800)
+    aZ = altar.properties.float(default=1000)
     aZ.doc = "the z semi-axis length"
 
     omegaX = altar.properties.float(default=0)
     omegaX.doc = "the CDM rotation about the x axis"
 
-    omegaY = altar.properties.float(default=-45)
+    omegaY = altar.properties.float(default=0)
     omegaY.doc = "the CDM rotation about the y axis"
 
     omegaZ = altar.properties.float(default=0)
     omegaZ.doc = "the CDM rotation about the z axis"
 
-    opening = altar.properties.float(default=1e2)
+    opening = altar.properties.float(default=50.0)
     opening.doc = "the tensile component of the Burgers vector of the dislocation"
 
     nu = altar.properties.float(default=.25)
@@ -106,7 +106,7 @@ def cdm(self):
         theta = π/4 # the azimuthal angle
         phi = π     # the polar angle
         # build the common projection vector
-        s = sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)
+        s = -sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta)
 
         # allocate a matrix to hold the projections
         los = altar.matrix(shape=(observations,3))
@@ -134,7 +134,7 @@ def cdm(self):
             # western stations
             if x < 0:
                 # come from a different data set
-                observation.oid = 1
+                observation.oid = 0
             # than
             else:
                 # eastern stations
@@ -158,7 +158,7 @@ def cdm(self):
         # go through the observations
         for idx, observation in enumerate(data):
             # set the covariance to a fraction of the "observed" displacement
-            correlation[idx,idx] = 1e-4
+            correlation[idx,idx] = 0.01**2
 
         # all done
         return data, correlation
@@ -171,9 +171,9 @@ def makeStations(self):
         # get some help
         import itertools
         # build a set of points on a grid
-        stations = itertools.product(range(-5,6), range(-5,6))
+        stations = itertools.product(range(-12,16), range(-12,16))
         # and return it
-        return tuple((x*1000, y*1000) for x,y in stations)
+        return tuple((x*500, y*500) for x,y in stations)
 
 
 # bootstrap

models/cdm/lib/libcdm/Source.cc

@@ -72,6 +72,9 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
     // clean up the resulting matrix
     gsl_matrix_set_zero(predicted);
 
+    // allocate storage for the displacement vectors; we reuse this for all samples
+    gsl_matrix * disp = gsl_matrix_alloc(_locations->size1, 3);
+
     // go through all the samples
     for (auto sample=0; sample<nSamples; ++sample) {
         // unpack the parameters
@@ -90,28 +93,62 @@ displacements(gsl_matrix_view * samples, gsl_matrix * predicted) const {
         auto omegaY = gsl_matrix_get(&samples->matrix, sample, _omegaYIdx);
         auto omegaZ = gsl_matrix_get(&samples->matrix, sample, _omegaZIdx);
 
-        // compute the displacements
-        cdm(sample, _locations, _los,
+    //     // compute the displacements
+    //     cdm(sample, _locations, _los,
+    //         xSrc, ySrc, dSrc,
+    //         aX, aY, aZ,
+    //         omegaX, omegaY, omegaZ,
+    //         openingSrc,
+    //         _nu,
+    //         predicted);
+
+    //     // apply the location specific projection to LOS vector and dataset shift
+    //     for (auto loc=0; loc<_locations->size1; ++loc) {
+    //         // get the current value
+    //         auto u = gsl_matrix_get(predicted, sample, loc);
+    //         // find the shift that corresponds to this observation
+    //         auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]);
+    //         // and apply it to the projected displacement
+    //         u -= shift;
+    //         // save
+    //         gsl_matrix_set(predicted, sample, loc, u);
+    //     }
+    // }
+        cdm(_locations,
             xSrc, ySrc, dSrc,
+            openingSrc,
             aX, aY, aZ,
             omegaX, omegaY, omegaZ,
-            openingSrc,
             _nu,
-            predicted);
-
+            disp);
+        
         // apply the location specific projection to LOS vector and dataset shift
         for (auto loc=0; loc<_locations->size1; ++loc) {
-            // get the current value
-            auto u = gsl_matrix_get(predicted, sample, loc);
+            // compute the components of the unit LOS vector
+            auto nx = gsl_matrix_get(_los, loc, 0);
+            auto ny = gsl_matrix_get(_los, loc, 1);
+            auto nz = gsl_matrix_get(_los, loc, 2);
+
+            // get the three components of the predicted displacement for this location
+            auto ux = gsl_matrix_get(disp, loc, 0);
+            auto uy = gsl_matrix_get(disp, loc, 1);
+            auto ud = gsl_matrix_get(disp, loc, 2);
+
+            // project
+            auto u = ux*nx + uy*ny + ud*nz;
             // find the shift that corresponds to this observation
             auto shift = gsl_matrix_get(&samples->matrix, sample, _offsetIdx+_oids[loc]);
             // and apply it to the projected displacement
             u -= shift;
+
             // save
             gsl_matrix_set(predicted, sample, loc, u);
         }
     }
 
+    // clean up
+    gsl_matrix_free(disp);
+
     // all done
     return;
 }