diff --git a/models/cdm/cdm/libcdm.py b/models/cdm/cdm/libcdm.py index 0ec13bbd..b3bc6b42 100644 --- a/models/cdm/cdm/libcdm.py +++ b/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: diff --git a/models/cdm/examples/cdm.pfg b/models/cdm/examples/cdm.pfg index 0b590f8d..53500f1d 100644 --- a/models/cdm/examples/cdm.pfg +++ b/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 diff --git a/models/cdm/examples/synthetic/cdm.py b/models/cdm/examples/synthetic/cdm.py index 7ad9a7e3..b282940f 100755 --- a/models/cdm/examples/synthetic/cdm.py +++ b/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 diff --git a/models/cdm/lib/libcdm/Source.cc b/models/cdm/lib/libcdm/Source.cc index 855b48b2..7786548e 100644 --- a/models/cdm/lib/libcdm/Source.cc +++ b/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; samplematrix, 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; } diff --git a/models/cdm/lib/libcdm/cdm.cc b/models/cdm/lib/libcdm/cdm.cc index 6bd26abe..a2fd9781 100644 --- a/models/cdm/lib/libcdm/cdm.cc +++ b/models/cdm/lib/libcdm/cdm.cc @@ -29,8 +29,12 @@ namespace altar { // local helpers static void - RDdispSurf(int sample, - const gsl_matrix * locations, const gsl_matrix * los, + // RDdispSurf(int sample, + // const gsl_matrix * locations, const gsl_matrix * los, + // const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, + // double opening, double nu, + // gsl_matrix * results); + RDdispSurf(const gsl_matrix * locations, const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, double opening, double nu, gsl_matrix * results); @@ -76,14 +80,21 @@ namespace altar { // definitions void altar::models::cdm:: -cdm(int sample, - const gsl_matrix * locations, const gsl_matrix * los, +cdm(const gsl_matrix * locations, double x, double y, double depth, + double opening, double aX, double aY, double aZ, double omegaX, double omegaY, double omegaZ, - double opening, double nu, gsl_matrix * predicted) +// cdm(int sample, +// const gsl_matrix * locations, const gsl_matrix * los, +// double x, double y, double depth, +// double aX, double aY, double aZ, +// double omegaX, double omegaY, double omegaZ, +// double opening, +// double nu, +// gsl_matrix * predicted) { // convert semi-axes to axes aX *= 2; @@ -91,12 +102,12 @@ cdm(int sample, aZ *= 2; // short circuit the trivial case - if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) { - // no displacements - gsl_matrix_set_zero(predicted); - // all done - return; - } + // if (std::abs(aX) < eps && std::abs(aY) < eps && std::abs(aZ) < eps) { + // // no displacements + // gsl_matrix_set_zero(predicted); + // // all done + // return; + // } // the axis specific coordinate matrices mat_t Rx = {1., 0., 0., @@ -141,22 +152,67 @@ cdm(int sample, Q1[2] > 0 || Q2[2] > 0 || Q3[2] > 0 || Q4[2] > 0 || R1[2] > 0 || R2[2] > 0 || R3[2] > 0 || R4[2] > 0) { // complain... - throw std::domain_error("the CDM must be below the surface"); + // throw std::domain_error("the CDM must be below the surface"); + printf("WARNING: The CDM must be below the surface \n"); } + // // dispatch the various cases + // if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) { + // RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); + // } else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) { + // RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); + // } else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) { + // RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // } else { + // RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); + // RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); + // RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // } + // dispatch the various cases - if (std::abs(aX) < eps && std::abs(aY) > eps && std::abs(aZ) > eps) { - RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); - } else if (std::abs(aX) > eps && std::abs(aY) < eps && std::abs(aZ) > eps) { - RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); - } else if (std::abs(aX) > eps && std::abs(aY) > eps && std::abs(aZ) < eps) { - RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + // short circuit the trivial case + if (std::abs(aX) == 0 && std::abs(aY) == 0 && std::abs(aZ) == 0) { + // no displacements + printf("I am if\n"); + gsl_matrix_set_zero(predicted); + } + else if (std::abs(aX) == 0 && std::abs(aY) > 0 && std::abs(aZ) > 0) { + printf("I am elseif1 \n"); + RDdispSurf(locations, P1, P2, P3, P4, opening, nu, predicted); + } else if (std::abs(aX) > 0 && std::abs(aY) == 0 && std::abs(aZ) > 0) { + printf("I am elseif2 \n"); + RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, predicted); + } else if (std::abs(aX) > 0 && std::abs(aY) > 0 && std::abs(aZ) == 0) { + printf("I am elseif3 \n"); + RDdispSurf(locations, R1, R2, R3, R4, opening, nu, predicted); } else { - RDdispSurf(sample, locations, los, P1, P2, P3, P4, opening, nu, predicted); - RDdispSurf(sample, locations, los, Q1, Q2, Q3, Q4, opening, nu, predicted); - RDdispSurf(sample, locations, los, R1, R2, R3, R4, opening, nu, predicted); + + // allocate room for the partial results + gsl_matrix * P = gsl_matrix_alloc(locations->size1, 3); + gsl_matrix * Q = gsl_matrix_alloc(locations->size1, 3); + gsl_matrix * R = gsl_matrix_alloc(locations->size1, 3); + + // compute + RDdispSurf(locations, P1, P2, P3, P4, opening, nu, P); + RDdispSurf(locations, Q1, Q2, Q3, Q4, opening, nu, Q); + RDdispSurf(locations, R1, R2, R3, R4, opening, nu, R); + + // assemble + for (auto loc=0; locsize1; ++loc) { + for (auto axis=0; axis<3; ++axis) { + // combine + auto result = + gsl_matrix_get(P, loc, axis) + + gsl_matrix_get(Q, loc, axis) + + gsl_matrix_get(R, loc, axis); + + // assign + gsl_matrix_set(predicted, loc, axis, result); + } + } } + // all done return; } @@ -164,11 +220,16 @@ cdm(int sample, // implementations static void altar::models::cdm:: -RDdispSurf(int sample, - const gsl_matrix * locations, const gsl_matrix * los, +// RDdispSurf(int sample, +// const gsl_matrix * locations, const gsl_matrix * los, +// const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, +// double opening, double nu, +// gsl_matrix * results) { +RDdispSurf(const gsl_matrix * locations, const vec_t & P1, const vec_t & P2, const vec_t & P3, const vec_t & P4, double opening, double nu, gsl_matrix * results) { + // cross auto V = cross(P2-P1, P4-P1); auto b = opening * V/norm(V); @@ -184,23 +245,30 @@ RDdispSurf(int sample, auto u3 = AngSetupFSC(x,y, b, P3, P4, nu); auto u4 = AngSetupFSC(x,y, b, P4, P1, nu); + // // assemble + // auto u = u1 + u2 + u3 + u4; + // // compute the unit LOS vector + // vec_t n = { gsl_matrix_get(los, loc, 0), + // gsl_matrix_get(los, loc, 1), + // gsl_matrix_get(los, loc, 2) }; + + // // project the displacement to the LOS + // auto uLOS = dot(u, n); + // // save by accumulating my contribution to the slot + // // N.B.: note the "+=": the general case call this function three times + // // get the current value + // auto current = gsl_matrix_get(results, sample, loc); + // // update + // current += uLOS; + // // save + // gsl_matrix_set(results, sample, loc, current); + // } + // assemble - auto u = u1 + u2 + u3 + u4; - // compute the unit LOS vector - vec_t n = { gsl_matrix_get(los, loc, 0), - gsl_matrix_get(los, loc, 1), - gsl_matrix_get(los, loc, 2) }; - - // project the displacement to the LOS - auto uLOS = dot(u, n); - // save by accumulating my contribution to the slot - // N.B.: note the "+=": the general case call this function three times - // get the current value - auto current = gsl_matrix_get(results, sample, loc); - // update - current += uLOS; - // save - gsl_matrix_set(results, sample, loc, current); + for (auto axis=0; axis<3; ++axis) { + auto u = u1[axis] + u2[axis] + u3[axis] + u4[axis]; + gsl_matrix_set(results, loc, axis, u); + } } // all done @@ -215,7 +283,8 @@ AngSetupFSC(double x, double y, double nu) { vec_t SideVec = PB - PA; vec_t eZ = {0, 0, 1}; - auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec)); + // auto beta = std::acos(dot(SideVec, eZ) / norm(SideVec)); + auto beta = std::acos(dot(-SideVec, eZ) / norm(SideVec)); if (std::abs(beta) < eps || std::abs(pi - beta) < eps) { return { 0,0,0 }; @@ -239,14 +308,20 @@ AngSetupFSC(double x, double y, vec_t vA, vB; // distinguish the two configurations - if (beta*adcsA[0] > 0) { + // if (beta*adcsA[0] > 0) { + if (beta*adcsA[0] >= 0) { // configuration I - vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]); - vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]); + // vA = AngDisDispSurf(adcsA, -pi+beta, b, nu, -PA[2]); + // vB = AngDisDispSurf(adcsB, -pi+beta, b, nu, -PB[2]); + vA = AngDisDispSurf(adcsA, -pi+beta, bADCS, nu, -PA[2]); + vB = AngDisDispSurf(adcsB, -pi+beta, bADCS, nu, -PB[2]); + } else { // configuration II - vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]); - vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]); + // vA = AngDisDispSurf(adcsA, beta, b, nu, -PB[2]); + // vB = AngDisDispSurf(adcsB, beta, b, nu, -PB[2]); + vA = AngDisDispSurf(adcsA, beta, bADCS, nu, -PA[2]); + vB = AngDisDispSurf(adcsB, beta, bADCS, nu, -PB[2]); } vec_t v = xform(transpose(A), vB - vA); diff --git a/models/cdm/lib/libcdm/cdm.h b/models/cdm/lib/libcdm/cdm.h index bfd6eb9a..906f665d 100644 --- a/models/cdm/lib/libcdm/cdm.h +++ b/models/cdm/lib/libcdm/cdm.h @@ -12,12 +12,20 @@ namespace altar { namespace models { namespace cdm { - void cdm(int sample, - const gsl_matrix * locations, const gsl_matrix * los, + // void cdm(int sample, + // const gsl_matrix * locations, const gsl_matrix * los, + // double X0, double Y0, double depth, + // double ax, double ay, double az, + // double omegaX, double omegaY, double omegaZ, + // double opening, + // double nu, + // gsl_matrix * predicted + // ); + void cdm(const gsl_matrix * locations, + double opening, double X0, double Y0, double depth, double ax, double ay, double az, double omegaX, double omegaY, double omegaZ, - double opening, double nu, gsl_matrix * predicted );