Skip to content

Commit

Permalink
Formatting cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
eheien committed Jul 8, 2014
1 parent f32417f commit 72e1f85
Show file tree
Hide file tree
Showing 8 changed files with 51 additions and 46 deletions.
3 changes: 2 additions & 1 deletion quakelib/src/QuakeLib.h
Original file line number Diff line number Diff line change
Expand Up @@ -163,9 +163,10 @@ namespace quakelib {
//! Set the maximum slip distance for this element
void set_max_slip(const double &new_max_slip) throw(std::invalid_argument) {
if (new_max_slip < 0) throw std::invalid_argument("quakelib::Element::set_max_slip");

_max_slip = new_max_slip;
}

//! Clear all variables for this element.
void clear(void) {
for (unsigned int i=0; i<3; ++i) {
Expand Down
40 changes: 21 additions & 19 deletions quakelib/src/QuakeLibIO.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -241,7 +241,7 @@ void quakelib::ModelElement::get_field_descs(std::vector<FieldDesc> &descs) {
field_desc.size = sizeof(float);
#endif
descs.push_back(field_desc);

field_desc.name = "max_slip";
field_desc.details = "Maximum slip distance for this element, in meters.";
#ifdef HDF5_FOUND
Expand Down Expand Up @@ -365,7 +365,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
double taper_full = 0;
double fault_area = 0;
cur_trace_point = element_start = conv.convert2xyz(trace.at(0).pos());

for (unsigned int i=1; i<trace.size(); ++i) {
next_trace_point = conv.convert2xyz(trace.at(i).pos());
t = 0;
Expand Down Expand Up @@ -395,19 +395,19 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac

// Use the t value in the middle of the element for interpolation
inter_t = t;
elem_depth = inter_t*trace.at(i).depth_along_dip()+(1.0-inter_t)*trace.at(i-1).depth_along_dip();
elem_depth = inter_t *trace.at(i).depth_along_dip()+(1.0-inter_t)*trace.at(i-1).depth_along_dip();
// Warn user if element depth is smaller than element size
num_vert_elems = floor(elem_depth/element_size);

if (num_vert_elems == 0) std::cerr << "WARNING: Depth is smaller than element size in trace segment "
<< i << ". Element size may be too big." << std::endl;

elem_slip_rate = conv.cm_per_yr2m_per_sec(inter_t*trace.at(i).slip_rate()+(1.0-inter_t)*trace.at(i-1).slip_rate());
elem_aseismic = inter_t*trace.at(i).aseismic()+(1.0-inter_t)*trace.at(i-1).aseismic();
elem_dip = conv.deg2rad(inter_t*trace.at(i).dip()+(1.0-inter_t)*trace.at(i-1).dip());
elem_rake = conv.deg2rad(inter_t*trace.at(i).rake()+(1.0-inter_t)*trace.at(i-1).rake());
elem_lame_mu = inter_t*trace.at(i).lame_mu()+(1.0-inter_t)*trace.at(i-1).lame_mu();
elem_lame_lambda = inter_t*trace.at(i).lame_lambda()+(1.0-inter_t)*trace.at(i-1).lame_lambda();
elem_slip_rate = conv.cm_per_yr2m_per_sec(inter_t *trace.at(i).slip_rate()+(1.0-inter_t)*trace.at(i-1).slip_rate());
elem_aseismic = inter_t *trace.at(i).aseismic()+(1.0-inter_t)*trace.at(i-1).aseismic();
elem_dip = conv.deg2rad(inter_t *trace.at(i).dip()+(1.0-inter_t)*trace.at(i-1).dip());
elem_rake = conv.deg2rad(inter_t *trace.at(i).rake()+(1.0-inter_t)*trace.at(i-1).rake());
elem_lame_mu = inter_t *trace.at(i).lame_mu()+(1.0-inter_t)*trace.at(i-1).lame_mu();
elem_lame_lambda = inter_t *trace.at(i).lame_lambda()+(1.0-inter_t)*trace.at(i-1).lame_lambda();

// Set up the vertical step to go along the dip
vert_step = element_step_vec.rotate_around_axis(Vec<3>(0,0,-1), M_PI/2);
Expand All @@ -418,20 +418,22 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
taper_t = 1;

double cur_dist = dist_along_trace+t*(next_trace_point-cur_trace_point).mag()-0.5*element_size;

if (taper_method == "taper_full" || taper_method == "taper_renorm") {
double x = cur_dist/total_trace_length;
double z = (float(ve)+0.5)/num_vert_elems;
taper_t *= 4*(x-x*x)*sqrt(1-z);
} else if (taper_method == "taper") {
double inside_dist = (total_trace_length/2.0 - fabs(total_trace_length/2.0-cur_dist));

if (inside_dist <= elem_depth) {
double x = inside_dist/elem_depth;
double z = (float(ve)+0.5)/num_vert_elems;
taper_t *= sqrt(x)*sqrt(1-z);
}
}

taper_flow += taper_t*elem_slip_rate*(element_size*element_size);
taper_flow += taper_t *elem_slip_rate*(element_size*element_size);
taper_full += elem_slip_rate*(element_size*element_size);

// Create the new vertices
Expand Down Expand Up @@ -462,7 +464,7 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
elem.set_rake(elem_rake);
elem.set_lame_mu(elem_lame_mu);
elem.set_lame_lambda(elem_lame_lambda);

fault_area += element_size*element_size;
}

Expand All @@ -485,11 +487,11 @@ void quakelib::ModelWorld::create_section(std::vector<unsigned int> &unused_trac
element(elem_ids[i]).set_slip_rate(renorm_factor*cur_slip_rate);
}
}

// Go through the created elements and assign maximum slip based on the total fault area
// From Table 2A in Wells Coppersmith 1994
double moment_magnitude = 4.07+0.98*log10(conv.sqm2sqkm(fault_area));

for (unsigned int i=0; i<elem_ids.size(); ++i) {
double max_slip = pow(10, (3.0/2.0)*(moment_magnitude+10.7))/(1e7*element(elem_ids[i]).lame_mu()*fault_area);
element(elem_ids[i]).set_max_slip(max_slip);
Expand Down Expand Up @@ -1080,8 +1082,8 @@ void quakelib::ModelWorld::write_section_hdf5(const hid_t &data_file) const {
ModelSection::get_field_descs(descs);
num_fields = descs.size();
num_sections = _sections.size();
field_names = new char*[num_fields];
field_details = new char*[num_fields];
field_names = new char *[num_fields];
field_details = new char *[num_fields];
field_offsets = new size_t[num_fields];
field_types = new hid_t[num_fields];
field_sizes = new size_t[num_fields];
Expand Down Expand Up @@ -1166,8 +1168,8 @@ void quakelib::ModelWorld::write_element_hdf5(const hid_t &data_file) const {
ModelElement::get_field_descs(descs);
num_fields = descs.size();
num_elements = _elements.size();
field_names = new char*[num_fields];
field_details = new char*[num_fields];
field_names = new char *[num_fields];
field_details = new char *[num_fields];
field_offsets = new size_t[num_fields];
field_types = new hid_t[num_fields];
field_sizes = new size_t[num_fields];
Expand Down Expand Up @@ -1252,8 +1254,8 @@ void quakelib::ModelWorld::write_vertex_hdf5(const hid_t &data_file) const {
ModelVertex::get_field_descs(descs);
num_fields = descs.size();
num_vertices = _vertices.size();
field_names = new char*[num_fields];
field_details = new char*[num_fields];
field_names = new char *[num_fields];
field_details = new char *[num_fields];
field_offsets = new size_t[num_fields];
field_types = new hid_t[num_fields];
field_sizes = new size_t[num_fields];
Expand Down
2 changes: 1 addition & 1 deletion quakelib/src/QuakeLibIO.h
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ namespace quakelib {
void set_lame_lambda(const float &lame_lambda) {
_data._lame_lambda = lame_lambda;
};

float max_slip(void) const {
return _data._max_slip;
};
Expand Down
2 changes: 1 addition & 1 deletion src/core/SimFramework.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
#include <iomanip>
#include <sstream>
#include <stdexcept>
#include <string.h>
#include <string.h>

#ifdef _OPENMP
#include <omp.h>
Expand Down
18 changes: 10 additions & 8 deletions src/core/VCSimulation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -548,19 +548,20 @@ void VCSimulation::distributeBlocks(const BlockIDSet &local_id_list, BlockIDProc
#ifdef DEBUG
stopTimer(dist_comm_timer);
#endif

// Create list of local block IDs
for (i=0,it=local_id_list.begin(); it!=local_id_list.end(); ++i,++it) local_block_ids[i] = *it;

// Count total, displacement of block IDs
int total_blocks = 0;
for (i=0;i<world_size;++i) {

for (i=0; i<world_size; ++i) {
proc_block_disps[i] = total_blocks;
total_blocks += proc_block_count[i];
}

BlockID *block_ids = new BlockID[total_blocks];

#ifdef DEBUG
startTimer(dist_comm_timer);
#endif
Expand All @@ -577,13 +578,14 @@ void VCSimulation::distributeBlocks(const BlockIDSet &local_id_list, BlockIDProc
#endif

n = 0;
for (p=0;p<world_size;++p) {
for (i=0;i<proc_block_count[p];++i) {

for (p=0; p<world_size; ++p) {
for (i=0; i<proc_block_count[p]; ++i) {
global_id_list.insert(std::make_pair(block_ids[n], p));
++n;
}
}

delete block_ids;
delete proc_block_count;
delete proc_block_disps;
Expand Down
2 changes: 1 addition & 1 deletion src/mesher.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,7 @@ void print_statistics(quakelib::ModelWorld &world, const std::string &file_name)
rake_vals.push_back(c.rad2deg(eit->rake()));
slip_rate_vals.push_back(c.m_per_sec2cm_per_yr(eit->slip_rate()));
}

if (num_elements == 0) {
rake_vals.push_back(std::numeric_limits<double>::quiet_NaN());
slip_rate_vals.push_back(std::numeric_limits<double>::quiet_NaN());
Expand Down
24 changes: 12 additions & 12 deletions src/simulation/RunEvent.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -134,15 +134,15 @@ void RunEvent::processBlocksSecondaryFailures(VCSimulation *sim, VCEventSweep &c

// Figure out how many failures there were over all processors
sim->distributeBlocks(local_id_list, global_id_list);

int num_local_failed = local_id_list.size();
int num_global_failed = global_id_list.size();

double *A = new double[num_local_failed*num_global_failed];
double *b = new double[num_local_failed];
double *x = new double[num_local_failed];

for (i=0,it=local_id_list.begin(); it!=local_id_list.end();++i,++it) {
for (i=0,it=local_id_list.begin(); it!=local_id_list.end(); ++i,++it) {
Block &blk = sim->getBlock(*it);

for (n=0,jt=global_id_list.begin(); jt!=global_id_list.end(); ++n,++jt) {
Expand All @@ -160,9 +160,9 @@ void RunEvent::processBlocksSecondaryFailures(VCSimulation *sim, VCEventSweep &c
double *fullA = new double[num_global_failed*num_global_failed];
double *fullb = new double[num_global_failed];
double *fullx = new double[num_global_failed];

// Fill in the A matrix and b vector from the various processes
for (i=0,n=0,jt=global_id_list.begin();jt!=global_id_list.end();++jt,++i) {
for (i=0,n=0,jt=global_id_list.begin(); jt!=global_id_list.end(); ++jt,++i) {
if (jt->second != sim->getNodeRank()) {
MPI_Recv(&(fullA[i*num_global_failed]), num_global_failed, MPI_DOUBLE, jt->second, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
MPI_Recv(&(fullb[i]), 1, MPI_DOUBLE, jt->second, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
Expand All @@ -172,31 +172,31 @@ void RunEvent::processBlocksSecondaryFailures(VCSimulation *sim, VCEventSweep &c
n++;
}
}

// Solve the global system on the root node
solve_it(num_global_failed, fullx, fullA, fullb);

// Send back the resulting values from x to each process
for (i=0,n=0,jt=global_id_list.begin();jt!=global_id_list.end();++jt,++i) {
for (i=0,n=0,jt=global_id_list.begin(); jt!=global_id_list.end(); ++jt,++i) {
if (jt->second != sim->getNodeRank()) {
MPI_Send(&(fullx[i]), 1, MPI_DOUBLE, jt->second, 0, MPI_COMM_WORLD);
} else {
memcpy(&(x[n]), &(fullx[i]), sizeof(double));
n++;
}
}

// Delete the memory arrays created
delete fullx;
delete fullb;
delete fullA;
} else {
for (i=0;i<num_local_failed;++i) {
for (i=0; i<num_local_failed; ++i) {
MPI_Send(&(A[i*num_global_failed]), num_global_failed, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
MPI_Send(&(b[i]), 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
}
for (i=0;i<num_local_failed;++i) {

for (i=0; i<num_local_failed; ++i) {
MPI_Recv(&(x[i]), 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}
}
Expand Down
6 changes: 3 additions & 3 deletions src/simulation/UpdateBlockStress.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -44,13 +44,13 @@ void UpdateBlockStress::init(SimFramework *_sim) {
// Set the stress drop based on the Greens function calculations
stress_drop = 0;
norm_velocity = sim->getBlock(gid).slip_rate();

for (nt=sim->begin(); nt!=sim->end(); ++nt) {
stress_drop += (nt->slip_rate()/norm_velocity)*sim->getGreenShear(gid, nt->getBlockID());
}

sim->getBlock(gid).setStressDrop(sim->getBlock(gid).max_slip()*stress_drop);

// noise in the range [1-a, 1+a)
sim->getBlock(gid).state.slipDeficit = 0;

Expand Down

0 comments on commit 72e1f85

Please sign in to comment.