From 4569a1fbb99d41ddd9027a65fdadc435627f7853 Mon Sep 17 00:00:00 2001 From: EstherHeck Date: Tue, 15 Mar 2022 12:03:32 +0100 Subject: [PATCH 1/2] add an extra erosional base level that prevents ghost nodes to lower too much --- include/aspect/mesh_deformation/fastscape.h | 43 +- source/mesh_deformation/fastscape.cc | 771 +++++++++++--------- 2 files changed, 436 insertions(+), 378 deletions(-) diff --git a/include/aspect/mesh_deformation/fastscape.h b/include/aspect/mesh_deformation/fastscape.h index 4c9dce78032..4fcf84b70bb 100644 --- a/include/aspect/mesh_deformation/fastscape.h +++ b/include/aspect/mesh_deformation/fastscape.h @@ -75,7 +75,7 @@ namespace aspect //void folder_output_(int *length, int *astep, const char *c); // Copy slopes, used in determined ghost node height for mass flux flow in from a boundary. void fastscape_copy_slope_(double *slopep); - + // View additional information from FastScape, not included in the .cc. void fastscape_view_(); void fastscape_debug_(); @@ -112,12 +112,12 @@ namespace aspect */ void parse_parameters (ParameterHandler &prm); - /** - * A function that fills the viscosity derivatives in the - * MaterialModelOutputs object that is handed over, if they exist. - * Does nothing otherwise. - */ - void set_ghost_nodes(double *h, double *vx, double *vy, double *vz, int nx, int ny) const; + /** + * A function that fills the viscosity derivatives in the + * MaterialModelOutputs object that is handed over, if they exist. + * Does nothing otherwise. + */ + void set_ghost_nodes(double *h, double *vx, double *vy, double *vz, int nx, int ny) const; private: // Number of FastScape steps per ASPECT timestep. @@ -213,7 +213,7 @@ namespace aspect double kdd; // Sediment transport coefficient. double kdsed; - + // Orographic parameters int mmax; int wb; @@ -222,6 +222,11 @@ namespace aspect double reduc_wb; bool stackoro; + // Parameters to set an extra erosional base level + // on the ghost nodes that differs from sea level. + bool use_extra_base_level; + double h_extra_base_level; + /** * Marine parameters */ @@ -259,18 +264,18 @@ namespace aspect double precision; - /** - * Interval between the generation of graphical output. This parameter - * is read from the input file and consequently is not part of the - * state that needs to be saved and restored. - */ - double output_interval; + /** + * Interval between the generation of graphical output. This parameter + * is read from the input file and consequently is not part of the + * state that needs to be saved and restored. + */ + double output_interval; - /** - * A time (in seconds) at which the last graphical output was supposed - * to be produced. Used to check for the next necessary output time. - */ - mutable double last_output_time; + /** + * A time (in seconds) at which the last graphical output was supposed + * to be produced. Used to check for the next necessary output time. + */ + mutable double last_output_time; }; } } diff --git a/source/mesh_deformation/fastscape.cc b/source/mesh_deformation/fastscape.cc index 1ad6a242493..aa1d0eaf6bf 100644 --- a/source/mesh_deformation/fastscape.cc +++ b/source/mesh_deformation/fastscape.cc @@ -65,31 +65,31 @@ namespace aspect // Initialize parameters for restarting fastscape restart = this->get_parameters().resume_computation; restart_step = 0; - - // Since we don't open these until we're on one processor, we need to check if the + + // Since we don't open these until we're on one processor, we need to check if the // restart files exist before hand. // TODO: This was quickly done and can likely be shortened/improved. - if(restart) - { - // Create variables for output directory and restart file - std::string dirname = this->get_output_directory(); - - std::ifstream in; - in.open(dirname + "fastscape_h_restart.txt"); - if (in.fail()) + if (restart) + { + // Create variables for output directory and restart file + std::string dirname = this->get_output_directory(); + + std::ifstream in; + in.open(dirname + "fastscape_h_restart.txt"); + if (in.fail()) AssertThrow(false,ExcMessage("Cannot open topography file to restart FastScape.")); - in.close(); - - in.open(dirname + "fastscape_b_restart.txt"); - if (in.fail()) + in.close(); + + in.open(dirname + "fastscape_b_restart.txt"); + if (in.fail()) AssertThrow(false,ExcMessage("Cannot open basement file to restart FastScape.")); - in.close(); - - in.open(dirname + "fastscape_steps_restart.txt"); - if (in.fail()) + in.close(); + + in.open(dirname + "fastscape_steps_restart.txt"); + if (in.fail()) AssertThrow(false,ExcMessage("Cannot open steps file to restart FastScape.")); - in.close(); - } + in.close(); + } // second is for maximum coordinates, first for minimum. for (unsigned int i=0; iget_time() >= last_output_time + output_interval || this->get_time()+a_dt >= end_time) - { - // Don't create a visualization file on a restart. - if(!restart) - make_vtk = 1; - - if (output_interval > 0) - { - // We need to find the last time output was supposed to be written. - // this is the last_output_time plus the largest positive multiple - // of output_intervals that passed since then. We need to handle the - // edge case where last_output_time+output_interval==current_time, - // we did an output and std::floor sadly rounds to zero. This is done - // by forcing std::floor to round 1.0-eps to 1.0. - const double magic = 1.0+2.0*std::numeric_limits::epsilon(); - last_output_time = last_output_time + std::floor((this->get_time()-last_output_time)/output_interval*magic) * output_interval/magic; + { + // Don't create a visualization file on a restart. + if (!restart) + make_vtk = 1; + + if (output_interval > 0) + { + // We need to find the last time output was supposed to be written. + // this is the last_output_time plus the largest positive multiple + // of output_intervals that passed since then. We need to handle the + // edge case where last_output_time+output_interval==current_time, + // we did an output and std::floor sadly rounds to zero. This is done + // by forcing std::floor to round 1.0-eps to 1.0. + const double magic = 1.0+2.0*std::numeric_limits::epsilon(); + last_output_time = last_output_time + std::floor((this->get_time()-last_output_time)/output_interval*magic) * output_interval/magic; + } } - } // Initialize kf and kd. for (int i=0; iget_time()/year_in_seconds; - double sediment_rain = sr_values[0]; - for (unsigned int j=0; jget_time()/year_in_seconds; + double sediment_rain = sr_values[0]; + for (unsigned int j=0; j sr_times[j]) - sediment_rain = sr_values[j+1]; + if (time > sr_times[j]) + sediment_rain = sr_values[j+1]; } /* @@ -466,7 +466,7 @@ namespace aspect const double h_seed = (std::rand()%( 2*noise_h+1 )) - noise_h; h[i] = h[i] + h_seed; } - + // Here we add the sediment rain (m/yr) as a flat increase in height. // This is done because adding it as an uplift rate would affect the basement. if (sediment_rain > 0 && use_marine) @@ -491,107 +491,107 @@ namespace aspect */ // I redid the indexing here, at some point I should double check that these still work without issue. if (use_ghost) - set_ghost_nodes(h.get(), vx.get(), vy.get(), vz.get(), nx, ny); + set_ghost_nodes(h.get(), vx.get(), vy.get(), vz.get(), nx, ny); //////// Section to apply orographic controls ///////// - // First for the wind barrier, we find the maximum height and index + // First for the wind barrier, we find the maximum height and index // along each line in the x and y direction. // If wind is east or west, we find maximum point for each ny row along x. std::vector> hmaxx(2, std::vector(ny, 0.0)); - if(wd == 0 || wd == 1) - { - for (int i=0; i hmaxx[0][i]) + if ( h[nx*i+j] > hmaxx[0][i]) { - hmaxx[0][i] = h[nx*i+j]; + hmaxx[0][i] = h[nx*i+j]; hmaxx[1][i] = j; } } } } - // If wind is north or south, we find maximum point for each nx row along y. - std::vector> hmaxy(2, std::vector(nx, 0.0)); - if(wd == 2 || wd == 3) - { + // If wind is north or south, we find maximum point for each nx row along y. + std::vector> hmaxy(2, std::vector(nx, 0.0)); + if (wd == 2 || wd == 3) + { for (int i=0; i hmaxy[0][i]) - { - hmaxy[0][i] = h[nx*j+i]; + if ( h[nx*j+i] > hmaxy[0][i]) + { + hmaxy[0][i] = h[nx*j+i]; hmaxy[1][i] = j; } } - } - } - + } + } + // Now we loop through all the points again and apply the reductions. for (int i=0; i wb) && (j < hmaxx[1][i]) || (h[nx*i+j] > mmax && !stackoro) ) - { - kf[nx*i+j] = kf[nx*i+j]*reduc_wb; - kd[nx*i+j] = kd[nx*i+j]*reduc_wb; - } - break; - } - case 1 : - { - if ( (hmaxx[0][i] > wb) && (j > hmaxx[1][i]) || (h[nx*i+j] > mmax && !stackoro) ) - { - kf[nx*i+j] = kf[nx*i+j]*reduc_wb; - kd[nx*i+j] = kd[nx*i+j]*reduc_wb; - } - break; - } - case 2 : + // Reduction from wind barrier. Apply a switch based off wind direction. + // Where 0 is wind going to the west, 1 the east, 2 the south, and 3 the north. + for (int j=0; j wb) && (i > hmaxy[1][j]) || (h[nx*i+j] > mmax && !stackoro) ) + switch (wd) { - kf[nx*i+j] = kf[nx*i+j]*reduc_wb; - kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + case 0 : + { + if ( (hmaxx[0][i] > wb) && (j < hmaxx[1][i]) || (h[nx*i+j] > mmax && !stackoro) ) + { + kf[nx*i+j] = kf[nx*i+j]*reduc_wb; + kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + } + break; + } + case 1 : + { + if ( (hmaxx[0][i] > wb) && (j > hmaxx[1][i]) || (h[nx*i+j] > mmax && !stackoro) ) + { + kf[nx*i+j] = kf[nx*i+j]*reduc_wb; + kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + } + break; + } + case 2 : + { + if ( (hmaxy[0][j] > wb) && (i > hmaxy[1][j]) || (h[nx*i+j] > mmax && !stackoro) ) + { + kf[nx*i+j] = kf[nx*i+j]*reduc_wb; + kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + } + break; + } + case 3 : + { + if ( (hmaxy[0][j] > wb) && (i < hmaxy[1][j]) || (h[nx*i+j] > mmax && !stackoro) ) + { + kf[nx*i+j] = kf[nx*i+j]*reduc_wb; + kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + } + break; + } + default : + AssertThrow(false, ExcMessage("This does not correspond with a wind direction.")); + break; } - break; - } - case 3 : - { - if ( (hmaxy[0][j] > wb) && (i < hmaxy[1][j]) || (h[nx*i+j] > mmax && !stackoro) ) + + // Apply elevation control. + if (h[nx*i+j] > mmax && stackoro) { - kf[nx*i+j] = kf[nx*i+j]*reduc_wb; - kd[nx*i+j] = kd[nx*i+j]*reduc_wb; + kf[nx*i+j] = kf[nx*i+j]*reduc_mmax; + kd[nx*i+j] = kd[nx*i+j]*reduc_mmax; } - break; - } - default : - AssertThrow(false, ExcMessage("This does not correspond with a wind direction.")); - break; - } - - // Apply elevation control. - if(h[nx*i+j] > mmax && stackoro) - { - kf[nx*i+j] = kf[nx*i+j]*reduc_mmax; - kd[nx*i+j] = kd[nx*i+j]*reduc_mmax; } - } - } - //////// End orographic section. Update to erosional parameters applied later. ///////// + } + //////// End orographic section. Update to erosional parameters applied later. ///////// // Get current fastscape timestep. int istep = 0; @@ -637,11 +637,11 @@ namespace aspect fastscape_set_dt_(&f_dt); // Set velocity components. - if(use_v) - { - fastscape_set_u_(vz.get()); - fastscape_set_v_(vx.get(), vy.get()); - } + if (use_v) + { + fastscape_set_u_(vz.get()); + fastscape_set_v_(vx.get(), vy.get()); + } // Set h to new values, and erosional parameters if there have been changes. fastscape_set_h_(h.get()); @@ -691,10 +691,10 @@ namespace aspect visualization_step = current_timestep; if (make_vtk) - { - this->get_pcout() << " Writing VTK..." << std::endl; - fastscape_named_vtk_(kd.get(), &vexp, &visualization_step, c, &length); - } + { + this->get_pcout() << " Writing VTK..." << std::endl; + fastscape_named_vtk_(kd.get(), &vexp, &visualization_step, c, &length); + } // If we've reached the end time, destroy fastscape. if (this->get_time()+a_dt >= end_time) @@ -835,240 +835,264 @@ namespace aspect template void FastScape::set_ghost_nodes(double *h, double *vx, double *vy, double *vz, int nx, int ny) const { - const int current_timestep = this->get_timestep_number (); - std::unique_ptr slopep (new double[array_size]()); - - /* - * Copy the slopes at each point, this will be used to set an H - * at the ghost nodes if a boundary mass flux is given. - */ - fastscape_copy_slope_(slopep.get()); - - /* - * Here we set the ghost nodes at the left and right boundaries. In most cases, - * this involves setting the node to the same values of v and h as the inward node. - * With the inward node being above or below for the bottom and top rows of ghost nodes, - * or to the left and right for the right and left columns of ghost nodes. - */ - for (int j=0; jget_timestep_number (); + std::unique_ptr slopep (new double[array_size]()); - /* - * Here we set the ghost nodes to the ones next to them, where for the left we - * add one to go to the node to the right, and for the right side - * we subtract one to go to the inner node to the left. - */ - vz[index_right] = vz[index_right-1]; - vz[index_left] = vz[index_left+1]; + /* + * Copy the slopes at each point, this will be used to set an H + * at the ghost nodes if a boundary mass flux is given. + */ + fastscape_copy_slope_(slopep.get()); - vy[index_right] = vy[index_right-1]; - vy[index_left] = vy[index_left+1]; + /* + * Here we set the ghost nodes at the left and right boundaries. In most cases, + * this involves setting the node to the same values of v and h as the inward node. + * With the inward node being above or below for the bottom and top rows of ghost nodes, + * or to the left and right for the right and left columns of ghost nodes. + */ + for (int j=0; j 0 && vx[index_left+1] >= 0) - { - side = index_right; - op_side = index_left; - jj = -1; - } - else if (vx[index_right-1] <= 0 && vx[index_left+1] < 0) - { - side = index_left; - op_side = index_right; - jj = 1; - } - else - continue; - - // Set right ghost node - h[index_right] = h[side+jj]; - vx[index_right] = vx[side+jj]; - vy[index_right] = vy[side+jj]; - vz[index_right] = vz[side+jj]; - - // Set left ghost node - h[index_left] = h[side+jj]; - vx[index_left] = vx[side+jj]; - vy[index_left] = vy[side+jj]; - vz[index_left] = vz[side+jj]; - - // Set opposing ASPECT boundary so it's periodic. - h[op_side-jj] = h[side+jj]; - vx[op_side-jj] = vx[side+jj]; - vz[op_side-jj] = vz[side+jj]; - vy[op_side-jj] = vy[side+jj]; - } - } + if (current_timestep == 1 || right_flux == 0) + { + slope = right_flux/kdd; + h[index_right] = h[index_right-1] + slope*2*dx; + if (use_extra_base_level) + { + if ((right==1) && (h[index_right] 0 && vx[index_left+1] >= 0) + { + side = index_right; + op_side = index_left; + jj = -1; + } + else if (vx[index_right-1] <= 0 && vx[index_left+1] < 0) + { + side = index_left; + op_side = index_right; + jj = 1; + } + else + continue; + + // Set right ghost node + h[index_right] = h[side+jj]; + vx[index_right] = vx[side+jj]; + vy[index_right] = vy[side+jj]; + vz[index_right] = vz[side+jj]; + + // Set left ghost node + h[index_left] = h[side+jj]; + vx[index_left] = vx[side+jj]; + vy[index_left] = vy[side+jj]; + vz[index_left] = vz[side+jj]; + + // Set opposing ASPECT boundary so it's periodic. + h[op_side-jj] = h[side+jj]; + vx[op_side-jj] = vx[side+jj]; + vz[op_side-jj] = vz[side+jj]; + vy[op_side-jj] = vy[side+jj]; + } + } - vx[index_bot] = vx[index_bot+nx]; - vx[index_top] = vx[index_top-nx]; + // Now do the same for the top and bottom ghost nodes. + for (int j=0; j 0 && vy[index_top-nx-1] >= 0) - { - side = index_top; - op_side = index_bot; - jj = -nx; - } - else if (vy[index_bot+nx-1] <= 0 && vy[index_top-nx-1] < 0) - { - side = index_bot; - op_side = index_top; - jj = nx; - } - else - continue; - - // Set top ghost node - h[index_top] = h[side+jj]; - vx[index_top] = vx[side+jj]; - vy[index_top] = vy[side+jj]; - vz[index_top] = vz[side+jj]; - - // Set bottom ghost node - h[index_bot] = h[side+jj]; - vx[index_bot] = vx[side+jj]; - vy[index_bot] = vy[side+jj]; - vz[index_bot] = vz[side+jj]; - - // Set opposing ASPECT boundary so it's periodic. - h[op_side-jj] = h[side+jj]; - vx[op_side-jj] = vx[side+jj]; - vz[op_side-jj] = vz[side+jj]; - vy[op_side-jj] = vy[side+jj]; - } - } + } + else + { + if (j == 0) + slope = top_flux/kdd - std::tan(slopep[index_top-nx+1]*numbers::PI/180.); + else if (j==(nx-1)) + slope = top_flux/kdd - std::tan(slopep[index_top-nx-1]*numbers::PI/180.); + else + slope = top_flux/kdd - std::tan(slopep[index_top-nx]*numbers::PI/180.); + + h[index_top] = h[index_top] + slope*2*dx; + } + + if (current_timestep == 1 || bottom_flux == 0) + { + slope = bottom_flux/kdd; + h[index_bot] = h[index_bot+nx] + slope*2*dx; + if (use_extra_base_level) + { + if ((bottom==1) && (h[index_bot] 0 && vy[index_top-nx-1] >= 0) + { + side = index_top; + op_side = index_bot; + jj = -nx; + } + else if (vy[index_bot+nx-1] <= 0 && vy[index_top-nx-1] < 0) + { + side = index_bot; + op_side = index_top; + jj = nx; + } + else + continue; + + // Set top ghost node + h[index_top] = h[side+jj]; + vx[index_top] = vx[side+jj]; + vy[index_top] = vy[side+jj]; + vz[index_top] = vz[side+jj]; + + // Set bottom ghost node + h[index_bot] = h[side+jj]; + vx[index_bot] = vx[side+jj]; + vy[index_bot] = vy[side+jj]; + vz[index_bot] = vz[side+jj]; + + // Set opposing ASPECT boundary so it's periodic. + h[op_side-jj] = h[side+jj]; + vx[op_side-jj] = vx[side+jj]; + vz[op_side-jj] = vz[side+jj]; + vy[op_side-jj] = vy[side+jj]; + } + } } // TODO: Give better explanations of variables and cite the fastscape documentation. @@ -1200,7 +1224,7 @@ namespace aspect prm.declare_entry("Sediment diffusivity", "-1", Patterns::Double(), "Diffusivity of sediment."); - prm.declare_entry("Orographic elevation control", "2000", + prm.declare_entry("Orographic elevation control", "2000", Patterns::Integer(), "Anything above this height has a change in erodibility."); prm.declare_entry("Orographic wind barrier height", "500", @@ -1213,11 +1237,32 @@ namespace aspect Patterns::Double(), "Amount to multiply kf and kd by past wind barrier."); prm.declare_entry ("Stack orographic controls", "false", - Patterns::Bool (), - "Whether or not to apply both controls to a point, or only a maximum of one set as the wind barrier."); + Patterns::Bool (), + "Whether or not to apply both controls to a point, or only a maximum of one set as the wind barrier."); prm.declare_entry ("Wind direction", "west", - Patterns::Selection("east|west|south|north"), - "This parameter assumes a wind direction, deciding which side is reduced from the wind barrier."); + Patterns::Selection("east|west|south|north"), + "This parameter assumes a wind direction, deciding which side is reduced from the wind barrier."); + prm.declare_entry ("Use an erosional base level differing from sea level", "false", + Patterns::Bool (), + "Whether or not to use an erosional base level that differs from sea level. Setting this parameter to " + "true will set all ghost nodes of fixed FastScape boundaries to the height you specify in " + "'set Erosional base level differing from sea level'. \nThis can make " + "sense for a continental model where the model surrounding topography is assumed above sea level, " + "e.g. highlands. If the sea level would be used as an erosional base level in this case, all topography " + "erodes away with lots of 'sediment volume' lost through the sides of the model. This is mostly " + "important, when there are mountains in the middle of the model, while it is less important when there " + "is lower relief in the middle of the model. \n" + "In the FastScape visualization files, setting the extra base level may show up as a strong " + "slope at the fixed boundaries of the model. However, in the ASPECT visualization files it will not " + "show up, as the ghost nodes only exist in FastScape."); + prm.declare_entry("Erosional base level differing from sea level", "0", + Patterns::Double(), + "When 'Use an erosional base level differing from sea level' is set to true, all ghost nodes of fixed " + "FastScape boundaries where no mass flux is specified by the user (FastScape boundary condition set to 1 " + "and 'Left/Right/Bottom/Top mass flux' set to 0) will be fixed to this elevation. The " + "reflecting boundaries (FastScape boundary condition set to 0) will not be affected, nor are the " + "boundaries where a mass flux is specified. \n" + "Units: m"); } prm.leave_subsection(); @@ -1287,12 +1332,12 @@ namespace aspect precision = prm.get_double("Precision"); noise_h = prm.get_integer("Initial noise magnitude"); sr_values = Utilities::string_to_double - (Utilities::split_string_list(prm.get ("Sediment rain"))); + (Utilities::split_string_list(prm.get ("Sediment rain"))); sr_times = Utilities::string_to_double - (Utilities::split_string_list(prm.get ("Sediment rain intervals"))); + (Utilities::split_string_list(prm.get ("Sediment rain intervals"))); - if (sr_values.size() != sr_times.size()+1) - AssertThrow(false, ExcMessage("Error: There must be one more sediment rain value than interval.")); + if (sr_values.size() != sr_times.size()+1) + AssertThrow(false, ExcMessage("Error: There must be one more sediment rain value than interval.")); prm.enter_subsection("Boundary conditions"); @@ -1331,18 +1376,26 @@ namespace aspect reduc_mmax = prm.get_double("Elevation factor"); reduc_wb = prm.get_double("Wind barrier factor"); stackoro = prm.get_bool("Stack orographic controls"); - - // Wind direction - if (prm.get ("Wind direction") == "west") - wd = 0; - else if (prm.get ("Wind direction") == "east") - wd = 1; - else if (prm.get ("Wind direction") == "north") - wd = 2; - else if (prm.get ("Wind direction") == "south") - wd = 3; - else - AssertThrow(false, ExcMessage("Not a valid wind direction.")); + + // Wind direction + if (prm.get ("Wind direction") == "west") + wd = 0; + else if (prm.get ("Wind direction") == "east") + wd = 1; + else if (prm.get ("Wind direction") == "north") + wd = 2; + else if (prm.get ("Wind direction") == "south") + wd = 3; + else + AssertThrow(false, ExcMessage("Not a valid wind direction.")); + + // set fixed ghost nodes to a base level for erosion that differs from sea level + use_extra_base_level = prm.get_bool("Use an erosional base level differing from sea level"); + if (use_extra_base_level) + AssertThrow(use_extra_base_level && use_ghost, ExcMessage( + "If you want to use an erosional base level differing from sea level, " + "you need to use ghost nodes.")); + h_extra_base_level = prm.get_double("Erosional base level differing from sea level"); } prm.leave_subsection(); From 0920ed85438bb66c1a6bbf7eec944d51848e93f9 Mon Sep 17 00:00:00 2001 From: EstherHeck Date: Tue, 15 Mar 2022 12:31:57 +0100 Subject: [PATCH 2/2] indent fastscape cc and h --- include/aspect/mesh_deformation/fastscape.h | 2 +- source/mesh_deformation/fastscape.cc | 46 ++++++++++----------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/include/aspect/mesh_deformation/fastscape.h b/include/aspect/mesh_deformation/fastscape.h index 4fcf84b70bb..41ec9a4710e 100644 --- a/include/aspect/mesh_deformation/fastscape.h +++ b/include/aspect/mesh_deformation/fastscape.h @@ -75,7 +75,7 @@ namespace aspect //void folder_output_(int *length, int *astep, const char *c); // Copy slopes, used in determined ghost node height for mass flux flow in from a boundary. void fastscape_copy_slope_(double *slopep); - + // View additional information from FastScape, not included in the .cc. void fastscape_view_(); void fastscape_debug_(); diff --git a/source/mesh_deformation/fastscape.cc b/source/mesh_deformation/fastscape.cc index aa1d0eaf6bf..41e67d8d6bc 100644 --- a/source/mesh_deformation/fastscape.cc +++ b/source/mesh_deformation/fastscape.cc @@ -885,14 +885,14 @@ namespace aspect */ slope = left_flux/kdd; h[index_left] = h[index_left+1] + slope*2*dx; - // This sets the ghost nodes of fixed boundaries to a user specified height which will - // be used as an erosional baselevel by fastscape and may limit erosional flux out - // of the model domain - if (use_extra_base_level) - { - if ((left==1) && (h[index_left]