From b42105db568d6ab39cec9ab2597e38c7b9780a3c Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 28 Jun 2022 09:59:51 +0100 Subject: [PATCH 1/2] Let gravprisms handle interfaces An interface grid is a surface that may cross zero. If so, the prisms we create should change sign across the zero threshold and for negative z we create a prism from x to 0 and for positive z we go from 0 to z. If a second surface or constant is given then we have a layer and no density changes sign. Documentation has ben updated to explain this as well adn a few comments were added elsewhere. I also added a test on the density contrast (-D cannot be 0, but variable density from grid could be zero at times and then we skip making a dummy prism). --- .../supplements/potential/gravprisms.rst | 23 ++++++- src/potential/gravprisms.c | 63 +++++++++++++------ 2 files changed, 64 insertions(+), 22 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index 45fd0f7a571..e44577f6f06 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -263,16 +263,33 @@ prism file and restrict calculations to the same crossing profile, i.e.:: gmt gravprisms -Ncrossing.txt -Mh @prisms.txt -Ff -Z7000 > faa_crossing.txt gmt plot faa_crossing.txt -R-30/30/0/350 -i0,3 -W1p -B -pdf faa_crossing +To build prisms using a variable density grid for an interface crossing the zero level +and obtain prisms with the negative density contrast if below zero and a positive +density contrast if above zero, then try:: -Note ----- + gmt gravprisms -TFlexure_surf.grd -C+wprisms_var.txt+q -DVariable_drho.grd + +Grids Straddling Zero Level +--------------------------- + +When creating prisms from grids via |-C|, a special case arises when a single surface (set via +|-L| or |-T|) straddles zero. This may happen if the surface reflects flexure beneath a +load, which has in a negative moat flanked by positive bulges. When such a *interface* grid +is detected we build prisms going from *z* to zero for negative *z* and from 0 to *z* for +positive *z*. As we flip below zero we also change the sign of the given density contrast. +You can override this behavior by specifying the opposite layer surface either by a constant +or another grid. E.g., if |-L| specifies the base surface you can eliminate prisms exceeding zero +via **-T**\ 0, and by interchanging the |-L| and |-T| arguments you can eliminate prisms below +zero. **Note**: When two surfaces are implied we keep the given density contrast as given. + +Note on Precision +----------------- The analytical expression for the geoid over a vertical prism (Nagy et al., 2000) is fairly involved and contains 48 terms. Due to various cancellations the end result is more unstable than the simpler expressions for gravity and VGG. Be aware that the result may have less significant digits that you may expect. - References ---------- diff --git a/src/potential/gravprisms.c b/src/potential/gravprisms.c index 130ef650610..4906456849f 100644 --- a/src/potential/gravprisms.c +++ b/src/potential/gravprisms.c @@ -203,8 +203,13 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT Ctrl->D.file = strdup (opt->arg); if (GMT_Get_FilePath (API, GMT_IS_GRID, GMT_IN, GMT_FILE_REMOTE, &(Ctrl->D.file))) n_errors++; } - else + else { /* Gave a constant density contrast */ Ctrl->D.rho = atof (opt->arg); + if (gmt_M_is_zero (Ctrl->D.rho)) { + GMT_Report (API, GMT_MSG_WARNING, "Option -D: Density contrast cannot be 0\n"); + n_errors++; + } + } break; case 'E': /* Set fixed prism dx, dy parameters instead of reading from prismfile */ n_errors += gmt_M_repeated_module_option (API, Ctrl->E.active); @@ -333,8 +338,8 @@ static int parse (struct GMT_CTRL *GMT, struct GRAVPRISMS_CTRL *Ctrl, struct GMT "Option -H: Cannot be used with -D.\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && !Ctrl->D.active && !Ctrl->H.active, "option -C: Need to select either -D or -H.\n"); - n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && Ctrl->C.dz == 0.0, - "Option -C: Requires +z when -H is selected\n"); + n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && Ctrl->C.dz <= 0.0, + "Option -C: Requires +z with a positive when -H is selected\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->C.active && Ctrl->H.active && !Ctrl->S.active, "Option -C: Requires -S when -H is set\n"); n_errors += gmt_M_check_condition (GMT, Ctrl->W.active && !Ctrl->C.active && !Ctrl->H.active, @@ -686,8 +691,9 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { GMT_Report (API, GMT_MSG_INFORMATION, "All z-values are assumed to be given in %s\n", uname[Ctrl->M.active[GRAVPRISMS_VER]]); if (Ctrl->C.active) { /* Need to create prisms from two surfaces first */ + bool flip_order = false; struct GMT_GRID *B = NULL, *T = NULL, *H = NULL, *Rho = NULL; - double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid, rs = 0.0, ws = 0.0; + double base = 0.0, top = 0.0, z1, z2, z_prev, z_next, z_mid, rs = 0.0, ws = 0.0, s = 1.0; size_t n_alloc = GMT_INITIAL_MEM_ROW_ALLOC; if (Ctrl->L.active) { /* Specified layer base */ @@ -719,53 +725,72 @@ EXTERN_MSC int GMT_gravprisms (void *V_API, int mode, void *args) { if ((Rho = GMT_Read_Data (API, GMT_IS_GRID, GMT_IS_FILE, GMT_IS_SURFACE, GMT_CONTAINER_AND_DATA, NULL, Ctrl->D.file, NULL)) == NULL) Return (API->error); } + + /* Some applications will give a interface that straddles z = 0. In that case we want prisms to + * go from a negative z to zero, and if z is positive, from 0 to a positive z. An example of such + * a surface is plate flexure deflections beneath a load, which will have positive bulges. In + * this case we also want to flip the sign of the given density contrast for negative z. */ + + if (Ctrl->L.active && !Ctrl->T.active && B) /* Just a base no top, does it straddle a zero level? */ + flip_order = (B->header->z_max > 0.0 && B->header->z_min < 0.0); + else if (Ctrl->T.active && !Ctrl->L.active && T) /* Just a top, no base, does it straddle a zero level? */ + flip_order = (T->header->z_max > 0.0 && T->header->z_min < 0.0); dx = 0.5 * H->header->inc[GMT_X]; dy = 0.5 * H->header->inc[GMT_Y]; rho = Ctrl->D.rho; /* Only initialized if -D is set */ - for (k = 0; k < 7; k++) + for (k = 0; k < 7; k++) /* Allocate memory for the 7 columns */ prism[k] = gmt_M_memory (GMT, NULL, n_alloc, double); - gmt_M_grd_loop (GMT, H, row, col, node) { - z2 = (T) ? T->data[node] : top; + gmt_M_grd_loop (GMT, H, row, col, node) { /* Visit every grid node */ + z2 = (T) ? T->data[node] : top; /* z-value at the top of the layer if (gmt_M_is_dnan (z2)) continue; /* Cannot work with NaNs - probably outside the feature */ - z1 = (B) ? B->data[node] : base; + z1 = (B) ? B->data[node] : base; /* z-value at the bottom of the layer */ if (gmt_M_is_dnan (z1)) continue; /* Cannot work with NaNs - probably outside the feature */ + if (Ctrl->D.file && gmt_M_is_zero (Rho->data[node])) continue; /* Need a nonzero density contrast to do any damage */ + if (flip_order) { /* E.g., when a flexure surface goes from moat to bulge */ + if (z2 < z1) { /* Must flip order and change density sign */ + gmt_M_double_swap (z1, z2); + s = -1.0; /* Flip density contrast too */ + } + else /* Good as is, reset density sign */ + s = 1.0; + } if (z2 <= z1) continue; /* No layer thickness detected */ z_prev = z1; /* We start exactly at z = z1 */ - do { /* Will at least be one prism */ - if (Ctrl->H.active) { /* Variable density means stacked prisms, most of thickness dz, at same point */ + do { /* There will at least be one prism */ + if (Ctrl->H.active) { /* Variable density model which means stacked prisms, most of thickness dz, at the same (x,y) */ /* First and last prisms in the stack are adjusted to have the height needed to match the surface boundaries */ z_next = (floor (z_prev / Ctrl->C.dz) + 1.0) * Ctrl->C.dz; /* Presumably next regular z-spacing */ if (z_next <= z_prev) z_next += Ctrl->C.dz; /* Can happen if z1 is a multiple of dz */ else if (z_next > z2) z_next = z2; /* At the top, clip to limit */ - z_mid = 0.5 * (z_prev + z_next); /* Middle of prism - used to look up density */ + z_mid = 0.5 * (z_prev + z_next); /* Middle of prism - used to look up density from model */ rho = gravprisms_mean_density (Ctrl, H->data[node], z_prev, z_next); } else { /* Constant density rho (set above via -D) or by Rho (via -W), just need a single prism per location */ z_next = z2; - if (Ctrl->D.file) - rho = Rho->data[node]; /* Update constant density for this prism */ + if (Ctrl->D.file) /* Update constant density for this prism */ + rho = Rho->data[node]; } if (n_prisms == n_alloc) { /* Need to allocate more memory for the prisms */ n_alloc <<= 1; /* Double it */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_alloc, double); } - /* Here we ensure prism is in meters */ + /* Here we ensure prism dimensions are given in meters */ prism[0][n_prisms] = scl_xy * (H->x[col] - dx); prism[1][n_prisms] = scl_xy * (H->x[col] + dx); prism[2][n_prisms] = scl_xy * (H->y[row] - dy); prism[3][n_prisms] = scl_xy * (H->y[row] + dy); prism[4][n_prisms] = scl_z * z_prev; prism[5][n_prisms] = scl_z * z_next; - prism[6][n_prisms] = rho; + prism[6][n_prisms] = s * rho; n_prisms++; - z_prev = z_next; /* The the top of this prism be the bottom of the next */ - } while (z_prev < z2); /* Until we run out of this stack */ + z_prev = z_next; /* The top of this prism will be the bottom of the next in the variable density case */ + } while (z_prev < z2); /* Until we run out of this stack of prisms */ if (Ctrl->W.active) { /* Get vertical average density and keep track of means */ - double dz = z2 - z1; + double dz = z2 - z1; /* Prism height */ Rho->data[node] = gravprisms_mean_density (Ctrl, H->data[node], z1, z2); rs += Rho->data[node] * dz; ws += dz; } } - /* Finalize allocation */ + /* Finalize memory allocation */ for (k = 0; k < 7; k++) prism[k] = gmt_M_memory (GMT, prism[k], n_prisms, double); GMT_Report (API, GMT_MSG_INFORMATION, "# of prisms constructed: %" PRIu64 "\n", n_prisms); From 3414eacb1fa0b3f58db446aa90b444a6c709bb35 Mon Sep 17 00:00:00 2001 From: Paul Wessel Date: Tue, 28 Jun 2022 10:05:15 +0100 Subject: [PATCH 2/2] Update gravprisms.rst --- doc/rst/source/supplements/potential/gravprisms.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/rst/source/supplements/potential/gravprisms.rst b/doc/rst/source/supplements/potential/gravprisms.rst index e44577f6f06..5fe51d36868 100644 --- a/doc/rst/source/supplements/potential/gravprisms.rst +++ b/doc/rst/source/supplements/potential/gravprisms.rst @@ -264,8 +264,8 @@ prism file and restrict calculations to the same crossing profile, i.e.:: gmt plot faa_crossing.txt -R-30/30/0/350 -i0,3 -W1p -B -pdf faa_crossing To build prisms using a variable density grid for an interface crossing the zero level -and obtain prisms with the negative density contrast if below zero and a positive -density contrast if above zero, then try:: +and obtain prisms with the negative of the given density contrast if below zero and the +positive density contrast if above zero, try:: gmt gravprisms -TFlexure_surf.grd -C+wprisms_var.txt+q -DVariable_drho.grd