egs: add support for MPI-based multi-component solves

1 files changed, 251 insertions, 0 deletions

diff --git a/components.c b/components.c
new file mode 100644
index 0000000..ef41653
--- /dev/null
+++ b/components.c
@@ -0,0 +1,251 @@
+/*
+ * Multi-component utility code
+ * Copyright 2019 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <errno.h>
+#include <stddef.h>
+#include <stdint.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include "common.h"
+#include "components.h"
+#include "mg2d_boundary.h"
+
+int mg2di_dg_copy(DomainGeometry **pdst, const DomainGeometry *src)
+{
+    DomainGeometry *dst;
+
+    dst = mg2di_dg_alloc(src->nb_components);
+    if (!dst)
+        return -ENOMEM;
+
+    memcpy(dst->components, src->components,
+           src->nb_components * sizeof(*src->components));
+    dst->domain_size[0] = src->domain_size[0];
+    dst->domain_size[1] = src->domain_size[1];
+
+    *pdst = dst;
+
+    return 0;
+}
+
+DomainGeometry *mg2di_dg_alloc(unsigned int nb_components)
+{
+    DomainGeometry *dg;
+
+    dg = calloc(1, sizeof(*dg));
+    if (!dg)
+        return NULL;
+
+    dg->components = calloc(nb_components, sizeof(*dg->components));
+    if (!dg->components)
+        goto fail;
+    dg->nb_components = nb_components;
+
+    return dg;
+fail:
+    mg2di_dg_free(&dg);
+    return NULL;
+}
+
+void mg2di_dg_free(DomainGeometry **pdg)
+{
+    DomainGeometry *dg = *pdg;
+
+    if (!dg)
+        return;
+
+    free(dg->components);
+
+    free(dg);
+    *pdg = NULL;
+}
+
+int mg2di_rect_intersect(Rect *dst, const Rect *src1, const Rect *src2)
+{
+    ptrdiff_t intersect_start0 = MAX(src1->start[0], src2->start[0]);
+    ptrdiff_t intersect_start1 = MAX(src1->start[1], src2->start[1]);
+    ptrdiff_t intersect_end0   = MIN(src1->start[0] + src1->size[0], src2->start[0] + src2->size[0]);
+    ptrdiff_t intersect_end1   = MIN(src1->start[1] + src1->size[1], src2->start[1] + src2->size[1]);
+
+    if (intersect_start0 < intersect_end0 && intersect_start1 < intersect_end1) {
+        dst->start[0] = intersect_start0;
+        dst->start[1] = intersect_start1;
+        dst->size[0]  = intersect_end0 - intersect_start0;
+        dst->size[1]  = intersect_end1 - intersect_start1;
+
+        return 1;
+    }
+
+    dst->size[0] = 0;
+    dst->size[1] = 0;
+
+    return 0;
+}
+
+/* merge adjacent/partially overlapping src into dst,
+ * if the result is rectangular */
+static int rect_merge(Rect *dst, const Rect *src)
+{
+    /* if dst is empty, copy src */
+    if (!dst->size[0] && !dst->size[1]) {
+        *dst = *src;
+        return 0;
+    }
+
+    /* if dst is fully inside src, copy src */
+    if (dst->start[0] >= src->start[0] &&
+        dst->start[0] + dst->size[0] <= src->start[0] + src->size[0] &&
+        dst->start[1] >= src->start[1] &&
+        dst->start[1] + dst->size[1] <= src->start[1] + src->size[1]) {
+        *dst = *src;
+        return 0;
+    }
+
+    /* if src is fully inside dst, do nothing */
+    if (src->start[0] >= dst->start[0] &&
+        src->start[0] + src->size[0] <= dst->start[0] + dst->size[0] &&
+        src->start[1] >= dst->start[1] &&
+        src->start[1] + src->size[1] <= dst->start[1] + dst->size[1])
+        return 0;
+
+    /* if src is adjacent or partially overlaps dst, merge them */
+    for (int dir = 0; dir < 2; dir++) {
+        if (dst->start[dir] != src->start[dir] ||
+            dst->size[dir]  != src->size[dir])
+            continue;
+
+        /* merge from above */
+        if (src->start[!dir] >= dst->start[!dir] &&
+            src->start[!dir] <= dst->start[!dir] + dst->size[!dir]) {
+            dst->size[!dir] = MAX(dst->size[!dir], src->start[!dir] + src->size[!dir] - dst->start[!dir]);
+            return 0;
+        }
+
+        /* merge from below */
+        if (dst->start[!dir] >= src->start[!dir] &&
+            dst->start[!dir] <= src->start[!dir] + src->size[!dir]) {
+
+            dst->size[!dir]  = MAX(src->size[!dir], dst->start[!dir] + dst->size[!dir] - src->start[!dir]);
+            dst->start[!dir] = src->start[!dir];
+            return 0;
+        }
+    }
+
+    /* src and dst do not match */
+    return -EINVAL;
+}
+
+static int calc_overlaps_component(Rect *overlaps, const DomainGeometry *dg,
+                                   unsigned int comp, int edge_width)
+{
+    const DomainComponent *c = &dg->components[comp];
+    int ret = 0;
+
+    Rect (*overlaps_tmp)[4];
+    int   *nb_overlaps_tmp;
+
+    overlaps_tmp    = calloc(dg->nb_components, sizeof(*overlaps_tmp));
+    nb_overlaps_tmp = calloc(dg->nb_components, sizeof(*nb_overlaps_tmp));
+    if (!overlaps_tmp || !nb_overlaps_tmp)
+        return -ENOMEM;
+
+    /* calculate the overlaps for each boundary layer */
+    for (int bnd_idx = 0; bnd_idx < 4; bnd_idx++) {
+        const int ci = mg2d_bnd_coord_idx(bnd_idx);
+        const int upper = mg2d_bnd_is_upper(bnd_idx);
+        Rect bnd;
+
+        bnd.start[ci] = c->interior.start[ci] + (upper ? c->interior.size[ci] : -edge_width);
+        bnd.size[ci]  = edge_width;
+
+        bnd.start[!ci] = c->interior.start[!ci] - edge_width;
+        bnd.size[!ci]  = c->interior.size[!ci] + 2 * edge_width;
+
+        for (unsigned int comp_idx = 0; comp_idx < dg->nb_components; comp_idx++) {
+            const DomainComponent *c1 = &dg->components[comp_idx];
+            Rect *dst = &overlaps_tmp[comp_idx][nb_overlaps_tmp[comp_idx]];
+
+            if (comp_idx == comp)
+                continue;
+
+            mg2di_rect_intersect(dst, &bnd, &c1->exterior);
+            if (dst->size[0] && dst->size[1])
+                nb_overlaps_tmp[comp_idx]++;
+        }
+    }
+
+    /* merge all the overlaps, we should get a single rectangle for each component */
+    for (unsigned int comp_idx = 0; comp_idx < dg->nb_components; comp_idx++) {
+        while (nb_overlaps_tmp[comp_idx] > 1) {
+            Rect *dst = &overlaps_tmp[comp_idx][0];
+            int merged = 0;
+
+            for (int i = 1; i < nb_overlaps_tmp[comp_idx]; i++) {
+                Rect *src = &overlaps_tmp[comp_idx][i];
+                ret = rect_merge(dst, src);
+                if (ret < 0)
+                    continue;
+
+                memmove(src, src + 1, nb_overlaps_tmp[comp_idx] - i - 1);
+                nb_overlaps_tmp[comp_idx]--;
+                merged = 1;
+                break;
+            }
+
+            if (!merged)
+                goto fail;
+        }
+        overlaps[comp_idx] = overlaps_tmp[comp_idx][0];
+    }
+
+fail:
+    free(overlaps_tmp);
+    free(nb_overlaps_tmp);
+
+    return ret;
+}
+
+int mg2di_dg_edge_overlaps(Rect *overlaps_recv, Rect *overlaps_send,
+                           const DomainGeometry *dg, unsigned int comp,
+                           unsigned int edge_width)
+{
+    Rect *overlaps = NULL;
+    int ret;
+
+    overlaps = calloc(SQR(dg->nb_components), sizeof(*overlaps));
+    if (!overlaps)
+        return -ENOMEM;
+
+    for (unsigned int i = 0; i < dg->nb_components; i++) {
+        ret = calc_overlaps_component(overlaps + i * dg->nb_components,
+                                      dg, i, edge_width);
+        if (ret < 0)
+            goto fail;
+    }
+
+    memcpy(overlaps_recv, overlaps + comp * dg->nb_components,
+           dg->nb_components * sizeof(*overlaps_recv));
+
+    for (unsigned int i = 0; i < dg->nb_components; i++)
+        overlaps_send[i] = overlaps[i * dg->nb_components + comp];
+
+fail:
+    free(overlaps);
+    return ret;
+}