(original) (raw)

struct ip { float *t_rp; /* time the reflected p-wave to reach nodes */ float *t_is; /* time the incident s-wave to reach nodes */ float *t_rs; /* time the reflected s-wave to reach nodes */ float c_s; /* speed of s-wave = sqrt(mv/E) */ float c_p; /* speed of p-wave */ float rho; /* speed of p-wave */ float cosG; /* cos(Gamma), Gamma denotes the plane */ /* incident wave comes along */ float sinG; /* sin(Gamma) */ float cosB; /* cos(Beta), angle between wave and the norm /* direction of the surface */ float sinB; /* sin(Beta) */ float cosA; /* cos(Alpha), reflected angle of the wave */ float sinA; /* sin(Alpha) */ float cAB; /* cos(Alpha+Beta) */ float kapa; /* c_s/c_p */ float d0; /* d0 = -x0*sinB + z0*cosB, x0,z0 are location of source of wave */ float f1; /* characteristic frequency of incident wave */ float ArBi; /* coefficient for computing free */ /* field motion */ float BrBi; /* coefficient for computing free */ /* field motion */ float dt; /* time step */ float duration; /* duration to be run, number of steps */ /* = duration/dt */ float dip; float strike; float rake; float T0; }; struct properties { float p_v; float s_v; float den; }; struct source { float dip, strike, rake; float x0[3]; }; struct source_node { int number; float x0[3]; float dist; }; struct triangle { float area; float cent[3]; int sign; }; static int ijk[4][3] = {{1,2,3},{2,3,0},{3,0,1},{0,1,2}}; #define PI 3.1415927 #define EPS 0.0001 #define COEF 0.017453293 #define DELTA 4 #define IDENTITY(x) x /* data types of output records */ #define OUT_CHAR 0 #define OUT_INT 1 #define OUT_FLOAT 2 struct source s; struct source_node sn; struct ip IP; void KMC_maker(); void prepare(); void iso3d(); void d_shape(); void inv_J(); float phi0(float); float phi1(float); float phi2(float); void cal_uvw(); void abc_damper(); void mv12x12(); void vv12x12(); float point2fault(float, float, float); float area(); void centroid(); void print_record_header(int, int); void print_file_header(); float dispx,dispy,dispz,velx,vely,velz; float ampl_h_dsp,ampl_h_vel,ampl_v_dsp,ampl_v_vel; float *max_h_dsp, *ang_h_dsp, *max_h_vel, *ang_h_vel; float *max_v_dsp, *max_v_vel, *energy_integral, *ampl_h_vel_prv; /* Define Archimedes vectors/matrices */ INODEVECTOR sourcenode; INODEVECTOR nodekind; NODEVECTOR nodekindf; IELEMVECTOR forcele; ELEMVECTOR elemkindf[3]; NODEVECTOR3 disp[3], vel, M, C; NODEVECTOR3 M23, C23, V23; MATRIX3 K; void main(int argc, char **argv) { int i, j, k, iter, timesteps; int NumNeighbor = 0; int NumPlus = 0; int NumMinus = 0; int globalNumNeighbor = 0; int globalNumPlus = 0; int disptplus, dispt, disptminus; int verticesonbnd; int cor[4], bv[3]; int Np_step; float c1[3], c2[3], c3[3], c4[3]; float time; float *Ke[12], Kee[12][12]; float Me[12], Ce[12], Mexv[12], Cexv[12], v[12]; float zeta, d_alpha, const_a, const_b; /* Rayleigh damping factors */ float x[3][4], xc[3]; float fr; /* Characteristic frequency for buffer zone */ float fault_u0; float uf[3]; struct properties prop; /* output related variables */ int FLOATS_PER_NODE = 7; /* output buffers */ int *ibuf; float *fbuf; int count; int surface_nodes; /* timing variables */ double starttime; double endtime; double totalstarttime; double totalendtime; double totaltime = 0.0; double totaltime1 = 0.0; double totaltime2 = 0.0; double packtime = 0.0; double readtime = 0.0; double printtime = 0.0; double inittime = 0.0; double sourcesetuptime = 0.0; double boundarytime = 0.0; double assembletime = 0.0; double mvtime = 0.0; double vvtime = 0.0; double misctime = 0.0; double mvflops; double mvmflops; double vvflops; double vvmflops; double localmflops; double globalmflops; double commbytes; double commthruput; /* NOTE: it has been assumed that there are 5 possible flag values for the node data, that is: 1 if the node is in the interior 4 if the node is on a x=const boundary surface 5 if the node is on a y=const boundary surface 6 if the node is on the bottom surface (z=z_lower) 3 if the node is on the surface (z=0) (but not along the edges) */ /*---------------------------------------------------------------*/ /* start timings */ TIME(totalstarttime); totalstarttime = totalstarttime - ARCHiotime; /* packfile loading time */ packtime = ARCHiotime; /* read data in from pack file */ TIME(ARCHstarttime); READNODEVECTOR(nodekindf); READELEMVECTOR(elemkindf[0]); READELEMVECTOR(elemkindf[1]); READELEMVECTOR(elemkindf[2]); TIME(ARCHendtime); ARCHiotime += (ARCHendtime - ARCHstarttime); readtime += (ARCHendtime - ARCHstarttime); TIME(ARCHstarttime); TIME(starttime); /* initializations */ /* dynamic parameters */ IP.dt = 0.0024; IP.duration = 40.0; IP.T0 = 0.6; Np_step = 30; timesteps = IP.duration / IP.dt + 1; /* source characteristics */ s.strike = 111.0; s.dip = 44.0; s.rake = 70.0; s.x0[0] = 32.264153; s.x0[1] = 23.814432; s.x0[2] = -11.25; IP.strike = s.strike * PI / 180.0; IP.dip = s.dip * PI / 180.0; IP.rake = s.rake * PI / 180.0; fault_u0 = 29.640788; /* prescribe slip motion */ uf[0] = uf[1] = uf[2] = 0.0; cal_uvw(&uf[0], &uf[1], &uf[2]); uf[0] *= fault_u0; uf[1] *= fault_u0; uf[2] *= fault_u0; /* damping constants */ zeta = 30.0; /* Damping factor FOR BUFFER ZONE in percentage(%) */ const_a = 0.00533333; const_b = 0.06666667; fr = 0.5; /* vector and array initialization */ M = 0.0; C = 0.0; M23 = 0.0; C23 = 0.0; V23 = 0.0; disp[0] = 0.0; disp[1] = 0.0; disp[2] = 0.0; nodekind = 0; K = 0.0; for (i = 0; i < 12; i++) Ke[i] = Kee[i]; disptplus = 0; dispt = 1; disptminus = 2; TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); inittime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* Print out case info */ if (ARCHself == 0) { fprintf(stderr, " CASE SUMMARY \n\n"); fprintf(stderr, " Fault Information\n"); fprintf(stderr, " Orientation: strike: %f\n",s.strike); fprintf(stderr, " dip: %f\n",s.dip); fprintf(stderr, " rake: %f\n",s.rake); fprintf(stderr, " Location: (%f, %f, %f), in Km\n", s.x0[0], s.x0[1], s.x0[2]); fprintf(stderr, " Dislocation: %f, cm\n", fault_u0); fprintf(stderr, " Ramp rise time: %f, sec\n", IP.T0); fprintf(stderr, " Time step: %f, sec\n", IP.dt); fprintf(stderr, " Duration: %f, sec\n", IP.duration); fflush(stderr); } TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* redefine nodekind to be 1 for all surface nodes */ /* this is temp work-around */ FORNODE(i) { nodekind[i] = (int) nodekindf[i]; if (nodekind[i] == 3) nodekind[i] = 1; } /* Count surface nodes */ surface_nodes = 0; for (i = 0; i < ARCHmine; i++) { GETCOORD3 (i,c1); if (c1[2] == 0.0) surface_nodes++; } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); inittime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); fprintf(stderr, "For subdomain %d there are %d surface nodes\n", ARCHself, surface_nodes); fflush(stderr); TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* various memory allocations */ if (surface_nodes != 0) { fbuf = (float *) malloc(surface_nodes * FLOATS_PER_NODE * sizeof(float)); if ( fbuf == (float *) NULL) { fprintf(stderr, "malloc failed for fbuf on processor = %d \n",ARCHself); fflush(stderr); exit(0); } ibuf = (int *) malloc(surface_nodes * sizeof(int)); if ( ibuf == (int *) NULL) { fprintf(stderr, "malloc failed for ibuf on processor = %d \n",ARCHself); fflush(stderr); exit(0); } max_h_dsp = (float *) malloc(surface_nodes * sizeof(float)); if ( max_h_dsp == (float *) NULL) { fprintf (stderr,"malloc failed for max_h_dsp on processor = %d \n",ARCHself); fflush(stderr); exit(0); } ang_h_dsp = (float *) malloc(surface_nodes * sizeof(float)); if ( ang_h_dsp == (float *) NULL) { fprintf (stderr,"malloc failed for ang_h_dsp on processor = %d \n",ARCHself); fflush(stderr); exit(0); } max_h_vel = (float *) malloc(surface_nodes * sizeof(float)); if ( max_h_vel == (float *) NULL) { fprintf (stderr,"malloc failed for max_h_vel on processor = %d \n",ARCHself); fflush(stderr); exit(0); } ang_h_vel = (float *) malloc(surface_nodes * sizeof(float)); if ( ang_h_vel == (float *) NULL) { fprintf (stderr,"malloc failed for ang_h_vel on processor = %d \n",ARCHself); fflush(stderr); exit(0); } max_v_dsp = (float *) malloc(surface_nodes * sizeof(float)); if ( max_v_dsp == (float *) NULL) { fprintf (stderr,"malloc failed for max_v_dsp on processor = %d \n",ARCHself); fflush(stderr); exit(0); } max_v_vel = (float *) malloc(surface_nodes * sizeof(float)); if ( max_v_vel == (float *) NULL) { fprintf (stderr,"malloc failed for max_v_vel on processor = %d \n",ARCHself); fflush(stderr); exit(0); } energy_integral = (float *) malloc(surface_nodes * sizeof(float)); if ( energy_integral == (float *) NULL) { fprintf (stderr,"malloc failed for energy_integral on processor = %d \n",ARCHself); fflush(stderr); exit(0); } ampl_h_vel_prv = (float *) malloc(surface_nodes * sizeof(float)); if ( ampl_h_vel_prv == (float *) NULL) { fprintf (stderr,"malloc failed for ampl_h_vel_prv on processor = %d \n",ARCHself); fflush(stderr); exit(0); } } else { fbuf = (float *) malloc (FLOATS_PER_NODE * sizeof(float)); ibuf = (int *) malloc (sizeof(int)); max_h_dsp = (float *) malloc(sizeof(float)); ang_h_dsp = (float *) malloc(sizeof(float)); max_h_vel = (float *) malloc(sizeof(float)); ang_h_vel = (float *) malloc(sizeof(float)); max_v_dsp = (float *) malloc(sizeof(float)); max_v_vel = (float *) malloc(sizeof(float)); energy_integral = (float *) malloc(sizeof(float)); ampl_h_vel_prv = (float *) malloc(sizeof(float)); } /* initialize output quantities */ for (i = 0; i < surface_nodes; i++) { max_h_dsp[i] = -1.0; max_h_vel[i] = -1.0; ang_h_dsp[i] = 0.0; ang_h_vel[i] = 0.0; max_v_dsp[i] = -1.0; max_v_vel[i] = -1.0; energy_integral[i] = 0.0; ampl_h_vel_prv[i] = 0.0; } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); inittime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* search for the node closest to the point source */ sourcenode = 1; forcele = 1; sn.number = 0; FORNODE(i) if (ARCHglobalnode[i] == 254216) { sn.number = i; sourcenode[i] = 10; } /* search for all the elements that contain the source node */ if (sn.number != 0) { FORELEM(i) { GETCORNERS(i, cor); if (cor[0] == sn.number || cor[1] == sn.number || cor[2] == sn.number || cor[3] == sn.number) { GETCOORD3(cor[0],c2); GETCOORD3(cor[1],c1); GETCOORD3(cor[2],c3); GETCOORD3(cor[3],c4); for(j = 0; j < 3; j++) { x[j][0] = c1[j]; x[j][1] = c2[j]; x[j][2] = c3[j]; x[j][3] = c4[j]; } centroid(x, xc); if (point2fault(xc[0], xc[1], xc[2]) >= 0) { NumPlus++; forcele[i] = 3; } else { NumMinus++; forcele[i] = 2; } NumNeighbor++; TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); sourcesetuptime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); fprintf(stderr, "ele %d at (%f %f %f) %d (subdomain %d)\n", i+1, xc[0], xc[1], xc[2], forcele[i],ARCHself); fprintf(stderr, "%d %d %d %d %d\n", i+1, cor[0], cor[1], cor[2], cor[3]); fflush(stderr); TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); } } } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); sourcesetuptime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); ISUMREDUCE(&NumNeighbor, &globalNumNeighbor, 1); ISUMREDUCE(&NumPlus, &globalNumPlus, 1); if (ARCHself == 0) { fprintf(stderr, "There are %d elements adjacent to the source\n", globalNumNeighbor); fprintf(stderr, " %d on the *plus* side of the mesh. \n", globalNumPlus); fflush(stderr); } TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* simulation */ FORELEM(i) { for (j = 0; j < 12; j++) { Me[j] = 0.0; Ce[j] = 0.0; v[j] = 0.0; for (k = 0; k < 12; k++) Kee[j][k] = 0.0; } GETCORNERS(i, cor); /* inhomogeneous material */ prop.p_v = elemkindf[0][i]; prop.s_v = elemkindf[1][i]; prop.den = elemkindf[2][i]; /* check if it is a boundary element */ TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); assembletime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); verticesonbnd = 0; for (j = 0; j < 4; j++) if (nodekind[cor[j]] != 1) bv[verticesonbnd++] = j; if (verticesonbnd == 3) { GETCOORD3(cor[bv[0]], c1); GETCOORD3(cor[bv[1]], c2); GETCOORD3(cor[bv[2]], c3); abc_damper(c1, c2, c3, bv, &prop, Ce); } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); boundarytime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); GETCOORD3(cor[0], c1); GETCOORD3(cor[1], c2); GETCOORD3(cor[2], c3); GETCOORD3(cor[3], c4); KMC_maker(c1, c2, c3, c4, &prop, Ke, Me, Ce); /* Here put damping into it with alpha only */ for(j = 0; j < 3; j++) { x[j][0] = c1[j]; x[j][1] = c2[j]; x[j][2] = c3[j]; x[j][3] = c4[j]; } centroid(x, xc); if (xc[2] < -11.5) d_alpha = 2.0 * zeta / 100.0 * (2.0 * PI * fr); else d_alpha = 4.0 * PI * const_a * 0.95 / (prop.s_v + const_b); for (j = 0; j < 12; j++) Ce[j] = Ce[j] + d_alpha * Me[j]; if (forcele[i] == 2 || forcele[i] == 3) { for (j = 0; j < 4; j++) { if (sourcenode[cor[j]] == 10) { v[3*j] = uf[0]; v[3*j+1] = uf[1]; v[3*j+2] = uf[2]; } else { v[3*j] = 0; v[3*j+1] = 0; v[3*j+2] = 0; } } vv12x12(Me, v, Mexv); vv12x12(Ce, v, Cexv); mv12x12(Kee, v); if (forcele[i] == 3) for (j = 0; j < 12; j++) { v[j] = -v[j]; Mexv[j] = -Mexv[j]; Cexv[j] = -Cexv[j]; } ASSEMBLEVECTOR3(v , i, V23); ASSEMBLEVECTOR3(Mexv, i, M23); ASSEMBLEVECTOR3(Cexv, i, C23); } ASSEMBLEVECTOR3(Me, i, M); ASSEMBLEVECTOR3(Ce, i, C); ASSEMBLEMATRIX3(Kee, i, K); } FULLASSEMBLE3(M); FULLASSEMBLE3(C); FULLASSEMBLE3(V23); FULLASSEMBLE3(M23); FULLASSEMBLE3(C23); TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); assembletime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); /* redefine nodekind to be 3 for all surface nodes */ count = 0; FORNODE(i) { if (i == ARCHmine) break; /* necessary to avoid overallocating ibuf,fbuf */ GETCOORD3 (i,c1); if (c1[2] == 0.0) { nodekind[i] = 3; ibuf[count] = ARCHglobalnode[i]; fbuf[2*count] = c1[0]; fbuf[2*count+1] = c1[1]; count++; /* on exit, count should be equal to surface_nodes */ } } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); inittime += (endtime - starttime); /* print the surface nodes ids and coordinates */ print_file_header(); TIME(ARCHstarttime); TIME(starttime); print_record_header(OUT_INT, surface_nodes); PRINTBUF(ibuf, surface_nodes*sizeof(int)); print_record_header(OUT_FLOAT, 2 * surface_nodes); PRINTBUF(fbuf, 2*surface_nodes*sizeof(float)); fbuf[0] = IP.dt; print_record_header(OUT_FLOAT, 1); PRINTBUF(fbuf, sizeof(float)); ibuf[0] = Np_step; ibuf[1] = timesteps; print_record_header(OUT_INT, 2); PRINTBUF(ibuf, 2*sizeof(int)); TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); /* time integration loop */ for (iter = 1; iter <= timesteps; iter++) { TIME(ARCHstarttime); TIME(starttime); MV3PRODUCT(K, disp[dispt], disp[disptplus]); TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); mvtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); time = iter * IP.dt; disp[disptplus] *= - IP.dt * IP.dt; disp[disptplus] += 2.0 * M * disp[dispt] - (M - IP.dt / 2.0 * C) * disp[disptminus] - IP.dt * IP.dt * (M23 * phi2(time) / 2.0 + C23 * phi1(time) / 2.0 + V23 * phi0(time) / 2.0); disp[disptplus] = disp[disptplus] / (M + IP.dt / 2.0 * C); vel = 0.5 / IP.dt * (disp[disptplus] - disp[disptminus]); TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); vvtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); j = 0; for (i = 0; i < ARCHmine; i++) { if (nodekind[i] == 3) { dispx = IDENTITY(disp[disptplus][i][0]); dispy = IDENTITY(disp[disptplus][i][1]); dispz = IDENTITY(disp[disptplus][i][2]); velx = IDENTITY(vel[i][0]); vely = IDENTITY(vel[i][1]); velz = IDENTITY(vel[i][2]); /* find maximum surface horizontal displacement; store amplitude and direction */ ampl_h_dsp = dispx * dispx + dispy * dispy; if (ampl_h_dsp > max_h_dsp[j]) { max_h_dsp[j] = ampl_h_dsp; ang_h_dsp[j] = atan2f(dispy,dispx); } /* find maximum surface horizontal velocity; store amplitude and direction */ ampl_h_vel = velx * velx + vely * vely; if (ampl_h_vel > max_h_vel[j]) { max_h_vel[j] = ampl_h_vel; ang_h_vel[j] = atan2f(vely,velx); } /* find maximum surface vertical displacement; store amplitude */ ampl_v_dsp = fabs(dispz); if (ampl_v_dsp > max_v_dsp[j]) max_v_dsp[j] = ampl_v_dsp; /* find maximum surface vertical velocity; store amplitude */ ampl_v_vel = fabs(velz); if (ampl_v_vel > max_v_vel[j]) max_v_vel[j] = ampl_v_vel; /* energy integral = int v^2 dt */ energy_integral[j] += 0.5 * (ampl_h_vel + ampl_h_vel_prv[j]) * IP.dt; ampl_h_vel_prv[j] = ampl_h_vel; j++; } } TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); misctime += (endtime - starttime); /* pack the output buffer and print */ TIME(ARCHstarttime); TIME(starttime); if (iter % Np_step == 0) { if (ARCHself == 0) { count = 1; fprintf(stderr,"Iteration = %d\n",iter); } else count = 0; ibuf[0] = iter; print_record_header(OUT_INT, count); PRINTBUF(ibuf, count*sizeof(int)); /* print response for all surface nodes */ count = 0; for (i = 0; i < ARCHmine; i++) { if (nodekind[i] == 3) { fbuf[count++] = IDENTITY(disp[disptplus][i][0]); fbuf[count++] = IDENTITY(disp[disptplus][i][1]); fbuf[count++] = IDENTITY(disp[disptplus][i][2]); fbuf[count++] = IDENTITY(vel[i][0]); fbuf[count++] = IDENTITY(vel[i][1]); fbuf[count++] = IDENTITY(vel[i][2]); } } print_record_header(OUT_FLOAT, count); PRINTBUF(fbuf, count*sizeof(float)); /* print response at the source node */ count = 0; for (i = 0; i < ARCHmine; i++) { if (sourcenode[i] == 10) { fbuf[count++] = IDENTITY(disp[disptplus][i][0]); fbuf[count++] = IDENTITY(disp[disptplus][i][1]); fbuf[count++] = IDENTITY(disp[disptplus][i][2]); fbuf[count++] = IDENTITY(vel[i][0]); fbuf[count++] = IDENTITY(vel[i][1]); fbuf[count++] = IDENTITY(vel[i][2]); } } print_record_header(OUT_FLOAT, count); PRINTBUF(fbuf, count*sizeof(float)); } TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(ARCHstarttime); TIME(starttime); i = disptminus; disptminus = dispt; dispt = disptplus; disptplus = i; TIME(ARCHendtime); TIME(endtime); ARCHcomptime += (ARCHendtime - ARCHstarttime); misctime += (endtime - starttime); } /* print maxima */ TIME(ARCHstarttime); TIME(starttime); count = 0; j = 0; for (i = 0; i < ARCHmine; i++) { if (nodekind[i] == 3) { fbuf[count++] = sqrtf(max_h_dsp[j]); fbuf[count++] = ang_h_dsp[j]; fbuf[count++] = sqrtf(max_h_vel[j]); fbuf[count++] = ang_h_vel[j]; fbuf[count++] = max_v_dsp[j]; fbuf[count++] = max_v_vel[j]; fbuf[count++] = energy_integral[j]; j++; } } print_record_header(OUT_FLOAT, count); PRINTBUF(fbuf, count*sizeof(float)); TIME(ARCHendtime); TIME(endtime); ARCHiotime += (ARCHendtime - ARCHstarttime); printtime += (endtime - starttime); TIME(totalendtime); /* various summaries of timings */ totaltime = (totalendtime - totalstarttime); totaltime1 = ARCHcomptime + ARCHcommtime + ARCHiotime; totaltime2 = packtime + readtime + printtime + inittime + sourcesetuptime + boundarytime + assembletime + mvtime + vvtime + misctime; /* flops and mflops per processor for matrix-vector operations */ mvflops = timesteps * 1.0 * /* 1 mvproduct/timestep */ (18 * /* 2 * dof^2 operations/entry */ (2*ARCHmatrixlen - ARCHnodes)); /* number of nonzero entries */ mvmflops = (mvflops/1000000.0)/mvtime; /* flops and mflops per processor for vector-vector operations */ vvflops = timesteps * 20.0 * /* twenty vv operations */ (ARCHnodes*3.0); /* three dof per operation */ vvmflops = (vvflops/1000000.0)/vvtime; /* total local mflops/processor (including all io and communication) */ localmflops = ((mvflops + vvflops)/1000000.0)/totaltime; /* aggregate global mflops (approximate because we assume */ /* equal work on each processor) */ globalmflops = (((mvflops + vvflops) * ARCHsubdomains) / 1000000.0) / totaltime; /* per processor communication during vector assembly */ commbytes = timesteps * 1.0 * /* one assembly ops per timestep */ ARCHcommlen * 6.0 * /* three items sent and received */ sizeof(float); /* bytes per item */ commthruput = (commbytes/(1<<20))/ARCHcommtime; if (ARCHself == 0) { fprintf(stderr, "%s: %d subdomains %d nodes %d elems %d timesteps\n", progname, ARCHsubdomains, ARCHglobalnodes, ARCHglobalelems, timesteps); fprintf(stderr, "\n"); fprintf(stderr, "comp : %12.2fs (%2d%%)\n", ARCHcomptime, (int)((ARCHcomptime/totaltime)*100.0)); fprintf(stderr, "comm : %12.2fs (%2d%%)\n", ARCHcommtime, (int)((ARCHcommtime/totaltime)*100.0)); fprintf(stderr, "io : %12.2fs (%2d%%)\n", ARCHiotime, (int)((ARCHiotime/totaltime)*100.0)); fprintf(stderr, "total : %12.2fs (%.0f min)[%.2f][%.2f]\n", totaltime, totaltime/60.0, totaltime1/totaltime, totaltime2/totaltime); fprintf(stderr, "total : %12.2fs\n",totaltime); fprintf(stderr, "total1 : %12.2fs\n",totaltime1); fprintf(stderr, "total2 : %12.2fs\n",totaltime2); fprintf(stderr, "pack file loading : %12.2fs (%2d%%)\n", packtime, (int)((packtime/totaltime)*100.0)); fprintf(stderr, "reading : %12.2fs (%2d%%)\n", readtime, (int)((readtime/totaltime)*100.0)); fprintf(stderr, "printing : %12.2fs (%2d%%)\n", printtime, (int)((printtime/totaltime)*100.0)); fprintf(stderr, "initializations : %12.2fs (%2d%%)\n", inittime, (int)((inittime/totaltime)*100.0)); fprintf(stderr, "source setup : %12.2fs (%2d%%)\n", sourcesetuptime, (int)((sourcesetuptime/totaltime)*100.0)); fprintf(stderr, "boundary setup : %12.2fs (%2d%%)\n", boundarytime, (int)((boundarytime/totaltime)*100.0)); fprintf(stderr, "assembling : %12.2fs (%2d%%)\n", assembletime, (int)((assembletime/totaltime)*100.0)); fprintf(stderr, "mv operations : %12.2fs (%2d%%)\n", mvtime, (int)((mvtime/totaltime)*100.0)); fprintf(stderr, "vv operations : %12.2fs (%2d%%)\n", vvtime, (int)((vvtime/totaltime)*100.0)); fprintf(stderr, "misc : %12.2fs (%2d%%)\n", misctime, (int)((misctime/totaltime)*100.0)); fprintf(stderr, "\n"); fprintf(stderr, "MFlops/sec/proc : %.1f\n", localmflops); fprintf(stderr, "MFlops/sec/proc (mv) : %.1f\n", mvmflops); fprintf(stderr, "MFlops/sec/proc (vv) : %.1f\n", vvmflops); fprintf(stderr, "Global MFlops/sec : %.1f\n", globalmflops); fprintf(stderr, "\n"); fprintf(stderr, "MBytes/sec/proc (as) : %.1f\n", commthruput); fprintf(stderr, "\n"); fprintf(stderr, "mvflops: %.0f mvtime: %.0f\n", mvflops, mvtime); fprintf(stderr, "vvflops: %.0f vvtime: %.0f\n", vvflops, vvtime); fprintf(stderr, "commbytes: %d\n", (int)commbytes); fflush(stderr); } } /* generates Ke[12][12], Me[12] and Ce[12] */ void KMC_maker(c1, c2, c3, c4, prop, Ke, Me, Ce) float *c1, *c2, *c3, *c4, **Ke, *Me, *Ce; struct properties *prop; { float jicob[3][5], ds[4][4], xyz[8][3], elas[4]; prepare(c1, c2, c3, c4, xyz, elas, prop); iso3d(xyz, elas, jicob, ds, Ke, Me); } /* returns material properties */ void Enu(e, v, prop) float *e, *v; struct properties *prop; { float x; x = prop->p_v / prop->s_v; x = x * x; *v = 0.5 * (x - 2) / (x - 1); *e = 2 * prop->den * prop->s_v * prop->s_v * (1 + *v); } /* returns corner coords and material properties */ void prepare(c1, c2, c3, c4, xyz, elas, prop) float *c1, *c2, *c3, *c4; float xyz[8][3]; float elas[]; struct properties *prop; { int i; float e, v; Enu(&e, &v, prop); elas[0] = e / (2.0 * (v + 1.0) * (1.0 - v * 2.0)); elas[1] = e * v / ((v + 1.0) * (1.0 - v * 2.0)); elas[2] = e / ((v + 1.0) * 2.0); elas[3] = prop->den; for (i = 0; i < 3; i++) { /* which defines corner coordinates */ xyz[0][i] = c1[i]; xyz[1][i] = c2[i]; xyz[2][i] = c3[i]; xyz[3][i] = c4[i]; } } /* generate element stiffness (Ke) and mass matrices (Me) for a 3d isoparameteric element (elasticity) stiffness matrix Ke and Me */ void iso3d(xyz, elas, jicob, ds, Ke, Me) float xyz[8][3]; float elas[]; float jicob[3][5]; float ds[4][4]; float **Ke; float *Me; { float det, tt, ts, c1, c2, c3; float volume; int i, j, m, n, row, column; d_shape(ds); /* compute derivertives of shapes */ for (i = 0; i < 3; ++i) { for (j = 0; j < 3; ++j) { det = 0.0; for (m = 0; m < 4; ++m) { det = det + ds[i][m]*xyz[m][j]; } jicob[j][i] = det; /* form Jacobian=dXj/dXIi */ } } inv_J(jicob, &det); /* find the J^-1 & its determinant */ for (m = 0; m < 4; ++m) { for (i = 0; i < 3; ++i) { jicob[i][3] = 0.0; for (j = 0; j < 3; ++j) { jicob[i][3] = jicob[i][3] + jicob[j][i]*ds[j][m]; } } for (i = 0; i < 3; ++i) { ds[i][m] = jicob[i][3]; /* d(shapes)/d(x,y,z) */ } } volume = det/6; /* volume of a tetrahedron in xyz coordinates */ /* printf("%f\n", volume); */ if (volume <= 0) { /* Something wrong w/ the numbering of ele vertices. */ printf("Warning: volume of the element = %f !\n", volume); /* exit(-1); */ } c1 = elas[0]*volume; /* volume in xi-coordinates: 1/6 */ c2 = elas[1]*volume; c3 = elas[2]*volume; row = -1; for (m = 0; m < 4; m++) { /* generate low-half matrix */ for (i = 0; i < 3; ++i) { ++row; column = -1; for (n = 0; n <= m; n++) { for (j = 0; j < 3; j++) { ++column; ts = ds[i][m]*ds[j][n]; if (i == j) { ts = ts*c1; tt = (ds[0][m]*ds[0][n] + ds[1][m]*ds[1][n] + ds[2][m]*ds[2][n])*c3; } else { if(m==n){ ts=ts*c1; tt=0; } else { ts = ts * c2; tt = ds[j][m] * ds[i][n] * c3; } } Ke[row][column] = Ke[row][column] + ts + tt; } } } } tt = elas[3] * volume / 4.0; for (i = 0; i < 12; i++) Me[i] = tt; for (i = 0; i < 12; ++i) for (j = 0; j <= i; j++) Ke[j][i] = Ke[i][j]; } /* calculate D(N_m)/D(XI_j) into ds[0..2][*] */ void d_shape(ds) float ds[4][4]; { ds[0][0] = -1; ds[1][0] = -1; ds[2][0] = -1; ds[0][1] = 1; ds[1][1] = 0; ds[2][1] = 0; ds[0][2] = 0; ds[1][2] = 1; ds[2][2] = 0; ds[0][3] = 0; ds[1][3] = 0; ds[2][3] = 1; } /* calculating the inverse and the determinant of the jacobian */ void inv_J(a, det) float a[3][5]; float *det; { float d1; float c[3][3]; int i,j; c[0][0] = a[1][1]*a[2][2] - a[2][1]*a[1][2]; c[0][1] = a[0][2]*a[2][1] - a[0][1]*a[2][2]; c[0][2] = a[0][1]*a[1][2] - a[0][2]*a[1][1]; c[1][0] = a[1][2]*a[2][0] - a[1][0]*a[2][2]; c[1][1] = a[0][0]*a[2][2] - a[0][2]*a[2][0]; c[1][2] = a[0][2]*a[1][0] - a[0][0]*a[1][2]; c[2][0] = a[1][0]*a[2][1] - a[1][1]*a[2][0]; c[2][1] = a[0][1]*a[2][0] - a[0][0]*a[2][1]; c[2][2] = a[0][0]*a[1][1] - a[0][1]*a[1][0]; *det = a[0][0]*c[0][0] + a[0][1]*c[1][0] + a[0][2]*c[2][0]; d1 = 1.0 / *det; for(i=0;i<3;i++) for(j=0;j<3;j++) a[i][j]= c[i][j]*d1; } /* calculate the area of a triangle given the coordinates of its three vertices (c1,c2,c3) */ float cal_area(c1, c2, c3) float *c1, *c2, *c3; { float a, b, c; float x2, y2, z2; float p; float area; x2 = (c1[0]-c2[0])*(c1[0]-c2[0]); y2 = (c1[1]-c2[1])*(c1[1]-c2[1]); z2 = (c1[2]-c2[2])*(c1[2]-c2[2]); a = sqrt(x2 + y2 + z2); x2 = (c3[0]-c2[0])*(c3[0]-c2[0]); y2 = (c3[1]-c2[1])*(c3[1]-c2[1]); z2 = (c3[2]-c2[2])*(c3[2]-c2[2]); b = sqrt(x2 + y2 + z2); x2 = (c1[0]-c3[0])*(c1[0]-c3[0]); y2 = (c1[1]-c3[1])*(c1[1]-c3[1]); z2 = (c1[2]-c3[2])*(c1[2]-c3[2]); c = sqrt(x2 + y2 + z2); p = (a + b + c) / 2; area = sqrt(p*(p-a)*(p-b)*(p-c)); if (area <= 0) { fprintf(stderr, "Error: negative area or zero area %f\n", p); exit(-1); } return area; } /* absorbing boundary element */ void abc_damper(c1, c2, c3, bv, prop, Ce) int bv[3]; float *c1, *c2, *c3; float Ce[12]; struct properties *prop; { int i, m0; float area; area = cal_area(c1, c2, c3); for (i = 0; i < 3; i++) { m0 = 3 * bv[i]; Ce[m0] = Ce[m0] + prop->s_v * prop->den * area / 3; Ce[m0+1] = Ce[m0+1] + prop->s_v * prop->den * area / 3; Ce[m0+2] = Ce[m0+2] + prop->p_v * prop->den * area / 3; } } /* ramp function */ float phi0(float t) { float value; if ( t <= IP.T0 ) { value = 0.5/PI*(2.0*PI*t/IP.T0-sinf(2.0*PI*t/IP.T0)); return value; } else return 1.0; } /* velocity of ramp function */ float phi1(float t) { float value; if (t <= IP.T0) { value = (1.0 - cosf(2.0 * PI * t / IP.T0)) / IP.T0; return value; } else return 0; } /* acceleration of ramp function */ float phi2(float t) { float value; if (t <= IP.T0) { value = 2 * PI / IP.T0 / IP.T0 * sinf(2 * PI * t / IP.T0); return value; } else return 0; } /* calculates the slip motion at the source node */ void cal_uvw(float *u, float *v, float *w) { *u = *v = *w = 0.0; *u = (cos(IP.rake)*sin(IP.strike)-sin(IP.rake)* cos(IP.strike)*cos(IP.dip)); *v = (cos(IP.rake)*cos(IP.strike)+sin(IP.rake)* sin(IP.strike)*cos(IP.dip)); *w = sin(IP.rake)*sin(IP.dip); } /* calculate the centroid of a tetrahedron */ void centroid(float x[3][4], float x0[3]) { int i; for (i = 0; i < 3; i++) x0[i] = (x[i][0] + x[i][1] + x[i][2] + x[i][3]) / 4; } /* calculate the distance to the fault from a given point (x,y,z) */ float point2fault(float x, float y, float z) { float nx, ny, nz; float d0; nx = cos(IP.strike) * sin(IP.dip); ny = - sin(IP.strike) * sin(IP.dip); nz = cos(IP.dip); d0 = - (nx * s.x0[0] + ny * s.x0[1] + nz * s.x0[2]); return (float) nx * x + ny * y + nz * z + d0; } /* matrix (12x12) times vector (12x1) product */ void mv12x12(float m[12][12], float v[12]) { int i, j; float u[12]; for (i = 0; i < 12; i++) { u[i] = 0; for (j = 0; j < 12; j++) u[i] += m[i][j] * v[j]; } for (i = 0; i < 12; i++) v[i] = u[i]; } /* vector (12x1) times vector (12x1) product */ void vv12x12 (float v1[12], float v2[12], float u[12]) { int i; for (i = 0; i < 12; i++) u[i] = v1[i] * v2[i]; } void print_record_header(int type, int items) { int rhdr[2], globalitems; ISUMREDUCE(&items, &globalitems, 1); if (ARCHself == 0) { rhdr[0] = type; rhdr[1] = globalitems; fwrite(rhdr, sizeof(int), 2, stdout); } } void print_file_header() { char buf[8]; buf[0] = 'A'; buf[1] = 'O'; buf[2] = '1'; if (ARCHself == 0) { fwrite(buf, sizeof(char), 8, stdout); } }