initial upload

fmm_2d_grid/fmm_2d_grid.cpp

@@ -0,0 +1,232 @@
+#include "iostream"
+#include "cmath"
+#include "vector"
+
+using namespace std;
+
+struct point_2d
+{
+	double x, y;
+};
+
+struct source : point_2d
+{
+	double rad;
+	double slow;
+};
+
+struct grid_point : point_2d
+{
+	int tag = 0; //0 = far away, 1 = close, 2 = active
+	double time = 1e+30;
+	double syn_time = 1e+30; //synthetic direct arrive time (only for error test)
+	double slow;
+	grid_point* neighbor[4] = {NULL,NULL,NULL,NULL}; //up down left right
+};
+typedef vector<grid_point> grid_point_array;
+typedef vector<grid_point*> ptr_grid_point_array;
+
+void update_heap(ptr_grid_point_array& a, int i, int n)
+{
+	int iMax = i,
+		iLeft = 2 * i + 1,
+		iRight = 2 * (i + 1);
+	if (iLeft < n && a[iMax]->time < a[iLeft]->time) {
+		iMax = iLeft;
+	}
+	if (iRight < n && a[iMax]->time < a[iRight]->time) {
+		iMax = iRight;
+	}
+	if (iMax != i) {
+		grid_point* tmp = a[iMax];
+		a[iMax] = a[i];
+		a[i] = tmp;
+		update_heap(a, iMax, n);
+	}
+	return;
+}
+
+void heap_sort(ptr_grid_point_array& a, int n)
+{
+	for (int i = (n - 1) / 2; i >= 0; i--) {
+		update_heap(a, i, n);
+	}
+	for (int i = n - 1; i > 0; --i) {
+		
+		grid_point* tmp = a[i];
+		a[i] = a[0];
+		a[0] = tmp;
+		update_heap(a, 0, i);
+	}
+	return;
+}
+
+double dis_point_2d(point_2d a, point_2d b)
+{
+	return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
+}
+
+void local_update(grid_point* point_ptr, double delta_h)
+{
+	//solve a quadratic equation to get the trial time
+	/*
+	double a = 0, b = 0, c = -1.0 * pow(delta_h,2) * pow(point_ptr->slow,2);
+	for (int i = 0; i < 4; i++)
+	{
+		if (point_ptr->neighbor[i] != NULL && point_ptr->neighbor[i]->tag == 2)
+		{
+			a += 1.0;
+			b += -2.0*point_ptr->neighbor[i]->time;
+			c += pow(point_ptr->neighbor[i]->time,2);
+		}
+	}
+	*/
+	// a second order upwind difference scheme
+	
+	double a = 0, b = 0, c = -4.0 * pow(delta_h,2) * pow(point_ptr->slow,2);
+	for (int i = 0; i < 4; i++)
+	{
+		if (point_ptr->neighbor[i] != NULL && point_ptr->neighbor[i]->tag == 2)
+		{
+			a += 4.0;
+			b += -8.0*point_ptr->neighbor[i]->time;
+			c += 4.0*pow(point_ptr->neighbor[i]->time,2);
+			if (point_ptr->neighbor[i]->neighbor[i] != NULL && point_ptr->neighbor[i]->neighbor[i]->tag == 2)
+			{
+				a += 5.0;
+				b += -16.0*point_ptr->neighbor[i]->time + 6.0*point_ptr->neighbor[i]->neighbor[i]->time;
+				c += 12.0*pow(point_ptr->neighbor[i]->time,2) - 8.0*point_ptr->neighbor[i]->time*point_ptr->neighbor[i]->neighbor[i]->time
+					+ pow(point_ptr->neighbor[i]->neighbor[i]->time,2);
+			}
+		}
+	}
+	
+	double delta = b*b - 4.0*a*c; //in a upwind scheme, delta is always bigger than zero
+	point_ptr->time = min(point_ptr->time, 0.5*(sqrt(delta) - b)/a);
+	return;
+}
+
+void abnormal_slowness(grid_point_array& grid_recall, double ab_slow,double xmin, double xmax, double ymin, double ymax)
+{
+	for (int i = 0; i < grid_recall.size(); i++)
+	{
+		if (grid_recall[i].x >= xmin && grid_recall[i].x <= xmax && grid_recall[i].y >= ymin && grid_recall[i].y <= ymax)
+		{
+			grid_recall[i].slow = ab_slow;
+		}
+	}
+	return;
+}
+
+int main(int argc, char const *argv[])
+{
+	//set grid parameters
+	int xnum = 201;
+	int ynum = 101;
+	double xmin = 0;
+	double ymin = 0;
+	double dh = 5;
+
+	//set source parameters
+	source init_source;
+	init_source.x = 50;
+	init_source.y = 250;
+	init_source.rad = 25;
+	init_source.slow = 1.0;
+
+	//initialize grid
+	grid_point_array grid_2d;
+	grid_2d.resize(xnum*ynum);
+	//down-left corner to up-right corner
+	for (int i = 0; i < ynum; i++)
+	{
+		for (int j = 0; j < xnum; j++)
+		{
+			grid_2d[i*xnum+j].x = xmin + dh*j;
+			grid_2d[i*xnum+j].y = ymin + dh*i;
+			grid_2d[i*xnum+j].slow = init_source.slow;
+
+			if (i <= ynum-2) grid_2d[i*xnum+j].neighbor[0] = &grid_2d[(i+1)*xnum+j]; //up
+			if (i >= 1) grid_2d[i*xnum+j].neighbor[1] = &grid_2d[(i-1)*xnum+j]; //down
+			if (j >= 1) grid_2d[i*xnum+j].neighbor[2] = &grid_2d[i*xnum+j-1]; //left
+			if (j <= xnum-2) grid_2d[i*xnum+j].neighbor[3] = &grid_2d[i*xnum+j+1]; //right
+
+			//calculate synthetic direct arrive time
+			grid_2d[i*xnum+j].syn_time = dis_point_2d(grid_2d[i*xnum+j], init_source) * init_source.slow;
+		}
+	}
+
+	//add abnormal slowness here
+	//abnormal_slowness(grid_2d,2.0,250,500,0,500);
+	//abnormal_slowness(grid_2d,4.0,800,1000,0,500);
+	//abnormal_slowness(grid_2d,1e+30,200,250,0,250);
+	//abnormal_slowness(grid_2d,1e+30,500,550,250,500);
+	//abnormal_slowness(grid_2d,1e+30,800,850,0,250);
+
+	ptr_grid_point_array close_node_ptr;
+	//initialize source nodes and close nodes;
+	for (int i = 0; i < xnum*ynum; i++)
+	{
+		if (dis_point_2d(grid_2d[i],init_source) <= init_source.rad)
+		{
+			grid_2d[i].tag = 2;
+			grid_2d[i].time = dis_point_2d(grid_2d[i],init_source) * init_source.slow;
+		}
+	}
+
+	for (int i = 0; i < xnum*ynum; i++)
+	{
+		if (grid_2d[i].tag == 2)
+		{
+			for (int j = 0; j < 4; j++)
+			{
+				if (grid_2d[i].neighbor[j] != NULL && grid_2d[i].neighbor[j]->tag == 0)
+				{
+					grid_2d[i].neighbor[j]->tag = 1;
+					close_node_ptr.push_back(grid_2d[i].neighbor[j]);
+				}
+			}
+		}
+	}
+
+	//calculate trial time for all close nodes
+	for (int i = 0; i < close_node_ptr.size(); i++)
+	{
+		local_update(close_node_ptr[i], dh);
+	}
+
+	//marching forward and updating the close nodes set
+	while (!close_node_ptr.empty())
+	{
+		// heap sort close nodes pointers to put the node first that has smallest time
+		heap_sort(close_node_ptr,close_node_ptr.size());
+
+		//change the first node's tag to 2 and update it's neighbor's time if it is not active
+		close_node_ptr[0]->tag = 2;
+		for (int i = 0; i < 4; i++)
+		{
+			if (close_node_ptr[0]->neighbor[i] != NULL && close_node_ptr[0]->neighbor[i]->tag == 0)
+			{
+				close_node_ptr[0]->neighbor[i]->tag = 1;
+				local_update(close_node_ptr[0]->neighbor[i], dh);
+				close_node_ptr.push_back(close_node_ptr[0]->neighbor[i]);
+			}
+			else if (close_node_ptr[0]->neighbor[i] != NULL && close_node_ptr[0]->neighbor[i]->tag == 1)
+			{
+				local_update(close_node_ptr[0]->neighbor[i], dh);
+			}
+		}
+		close_node_ptr.erase(close_node_ptr.begin());
+	}
+
+	for (int i = 0; i < xnum*ynum; i++)
+	{
+		if (grid_2d[i].time > 1e+10)
+		{
+			cout << grid_2d[i].x << " " << grid_2d[i].y << " " << grid_2d[i].tag << " " << grid_2d[i].syn_time << " nan " << grid_2d[i].time - grid_2d[i].syn_time << endl;
+		}
+		else cout << grid_2d[i].x << " " << grid_2d[i].y << " " << grid_2d[i].tag << " " << grid_2d[i].syn_time << " " << grid_2d[i].time << " " << grid_2d[i].time - grid_2d[i].syn_time << endl;
+	}
+	
+	return 0;
+}

fmm_2d_grid/run.sh

@@ -0,0 +1,9 @@
+#!/bin/bash
+
+g++-9 fmm_2d_grid.cpp -o fmm_2d_grid -O2
+
+./fmm_2d_grid > out.txt
+
+xyz2grd out.txt -Gout.nc -R0/1000/0/500 -I5/5 -i0,1,5
+
+ncview out.nc