parallel second order FSM

Author: Dengda98

docs/source/examples/01_uniform.ipynb

@@ -15,8 +15,10 @@ Jupyter 使用示例
    :maxdepth: 2
 
    examples/01_uniform
+   examples/01_uniform_spherical
    examples/02_random_model
    examples/02_cplxmodel_raytracing
    examples/03_spherical
+   examples/03_spherical2
    examples/04_distance
 

docs/source/examples/01_uniform_spherical.ipynb

@@ -35,8 +35,9 @@ LIB_NAME = $(LIB_DIR)$(FILESEP)libfmm
 LIB_EXT = .so
 
 CC = gcc
-CFLAGS =  -O3 -g -ffast-math -march=native -mtune=native -fPIC -Wimplicit-function-declaration -I$(INC_DIR) -lm   # -fsanitize=address -lasan
-LDFLAGS = -shared -lm # -lasan
+FOMPFLAGS = -fopenmp 
+CFLAGS = $(FOMPFLAGS) -O3 -g -ffast-math -march=native -mtune=native -fPIC -Wimplicit-function-declaration -I$(INC_DIR) -lm   # -fsanitize=address -lasan
+LDFLAGS = $(FOMPFLAGS) -shared -lm # -lasan
 
 SRCS = $(wildcard $(SRC_DIR)/*.c)
 INCS = $(wildcard $(INC_DIR)/*.h)

docs/source/examples/02_cplxmodel_raytracing.ipynb

@@ -7,9 +7,79 @@
 #pragma once
 
 #include "const.h"
+#include "heapsort.h"
+
+
+/**
+ * 定义OpenMP多线程数
+ * 
+ * @param     num_threads     (in)线程数
+ */
+void set_fsm_num_threads(int num_threads);
+
+
+/**
+ * 使用Fast Sweeping Method计算全局走时场
+ * 
+ * @param     rs     (in)维度1坐标数组
+ * @param     nr     (in)rs长度
+ * @param     ts     (in)维度2坐标数组
+ * @param     nt     (in)ts长度
+ * @param     ps     (in)维度2坐标数组
+ * @param     np     (in)ps长度
+ * @param     rr     (in)源点维度1坐标
+ * @param     tt     (in)源点维度2坐标
+ * @param     pp     (in)源点维度3坐标
+ * @param     maxodr (in)使用的最大差分阶数
+ * @param     Slw    (in)展平的三维慢度场
+ * @param     TT     (inout)展平的三维走时场，如果初始值有非零值，会被直接加入堆中，此时源点不再使用
+ * @param     sphcoord  (in)是否使用球坐标
+ * @param     rfgfac    (in)对于源点附近的格点间加密倍数，>1
+ * @param     rfgn      (in)对于源点附近的格点间加密处理的辐射半径，>=1
+ * @param     printbar  (in)是否打印进度条
+ * @param     eps        (in)Sweep后的最大更新量达到收敛条件
+ * @param     maxLoops   (in)Fast Sweeping Method整体迭代次数（对于3D模型，向8个方向各Sweep一次为迭代一次）  
+ * @param     isparallel (in)是否使用并行FSM
+ * 
+ * 
+ */
+int FastSweeping(
+    const double *rs, int nr, 
+    const double *ts, int nt, 
+    const double *ps, int np,
+    double rr,  double tt, double pp,
+    int maxodr,  const MYREAL *Slw, 
+    MYREAL *TT, bool sphcoord, 
+    int rfgfac, int rfgn, bool printbar, 
+    double eps, int maxLoops, bool isparallel);
+
+
+/**
+ * 在有初始走时的情况下使用Fast Marching Method计算全局走时场
+ * 
+ * @param     rs     (in)维度1坐标数组
+ * @param     nr     (in)rs长度
+ * @param     ts     (in)维度2坐标数组
+ * @param     nt     (in)ts长度
+ * @param     ps     (in)维度2坐标数组
+ * @param     np     (in)ps长度
+ * @param     maxodr (in)使用的最大差分阶数
+ * @param     Slw    (in)展平的三维慢度场
+ * @param     TT     (inout)展平的三维走时场
+ * @param     FMM_stat  (out)记录每个节点的状态(alive, close, far)
+ * @param     sphcoord  (in)是否使用球坐标
+ * @param     printbar  (in)是否打印进度条
+ * @param     eps        (in)Sweep后的最大更新量达到收敛条件
+ * @param     maxLoops   (in)Fast Sweeping Method整体迭代次数（对于3D模型，向8个方向各Sweep一次为迭代一次）  
+ * @param     isparallel (in)是否使用并行FSM
+ * 
+ * @return    nsweep, 
+ */
+int FastSweeping_with_initial(
+    const double *rs, int nr, 
+    const double *ts, int nt, 
+    const double *ps, int np,
+    int maxodr,  const MYREAL *Slw, MYREAL *TT, 
+    char *FMM_stat, bool sphcoord, bool printbar, 
+    double eps, int maxLoops, bool isparallel);
 
-typedef struct _FSW_DATA {
-    // int idx;
-    MYREAL t;
-    double h;
-} FSW_DATA;

docs/source/examples/02_random_model.ipynb

@@ -137,8 +137,8 @@ HEAP_DATA * FastMarching_with_initial(
     double sin_ts[nt];
     if(sphcoord){
         for(int it=0; it<nt; ++it){
-            sin_ts[it] = sin(ts[it]);
-            if(fabs(sin_ts[it]) < 1e-12) sin_ts[it] += 1e-12;
+            sin_ts[it] = fabs(sin(ts[it]));
+            if(sin_ts[it] < 1e-12) sin_ts[it] += 1e-12;
         }
     }
 
@@ -592,7 +592,7 @@ HEAP_DATA * init_source_TT_refinegrid(
             TT[jdx] = rfg_TT[rfg_jdx];
             // printf("rfg_TT = %f, jr, jt, jp %d, %d, %d\n", rfg_TT[rfg_jdx], jr, jt, jp);
             FMM_stat[jdx] = FMM_ALV;
-            (*pNdots)--;
+            if(pNdots!=NULL) (*pNdots)--;
         }
         if(jdx1>=0){
             FMM_data = HeapPush(FMM_data, psize, pcap, jdx1, NroIdx, TT);
@@ -828,16 +828,6 @@ MYREAL FMM_raytracing(
 
 
     int ntp = nt*np;
-
-    double sin_ts[nt];
-    if(sphcoord){
-        for(int it=0; it<nt; ++it){
-            sin_ts[it] = sin(ts[it]);
-            if(fabs(sin_ts[it]) < 1e-6) sin_ts[it] += 1e-6;
-        }
-    } else {
-        for(int it=0; it<nt; sin_ts[it]=0.0, ++it);
-    }
 
     int idot = 0;