diff --git a/CMakeLists.txt b/CMakeLists.txt index 815006d2947..72b0eacafc1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -262,6 +262,8 @@ else() message(STATUS "Using ASPECT_WITH_NETCDF = 'OFF'") endif() +# openLEM +include_directories("${CMAKE_SOURCE_DIR}/contrib/openlem/") # WorldBuilder set(ASPECT_WITH_WORLD_BUILDER ON CACHE BOOL "Whether to enable compiling aspect with the Geodynamic World Builder.") diff --git a/contrib/openlem/openlem.cpp b/contrib/openlem/openlem.cpp new file mode 100644 index 00000000000..6058d37fa61 --- /dev/null +++ b/contrib/openlem/openlem.cpp @@ -0,0 +1,3310 @@ +// OpenLEM Version 44 experimental 2025-07-10 +// +// Copyright (C) 2012-2025 Stefan Hergarten +// +// This program is free software; you can redistribute it and/or modify it +// under the terms of the GNU General Public License as published by the +// Free Software Foundation; either version 3 of the License, or (at your +// option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +// for more details. To view a copy of this license, visit +// https://www.gnu.org/licenses/gpl-3.0.txt or write to the Free Software +// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. +// +// For any questions concerning the codes and bug reports, please contact the +// developer at stefan.hergarten@geologie.uni-freiburg.de. +// +// Changelog +// +// Version 44 2025-07-10 +// Reimplemented the method fillLakes for better efficiency. +// Fixed problem with flow directions if boundaries are changed during the +// simulation. +// Changed status of flooded areas (water level > topography) from hillslope +// to channel for option CHANNEL. +// +// Version 43 2023-08-22 +// Introduced faster computation for neighbors. +// Changed some default parameter values. +// Renamed method computeDischarge to computeFluxes. +// Removed option INTOLAKES. +// Introduced new concept for computing the flow pattern in lakes (option +// LAKES without LAKEFLOWDIR, highly experimental). +// Introduced option REDREC for reducing the level of recursion (highly +// experimental) +// +// Version 42 2023-06-29 +// Fixed problem with stream power in lakes (option LAKES). +// Introduced option BORROW for treating the situation that the erosion +// of deposits exceeds the available thickness of the deposit layer +// (in combination with the option LAYERS). This option removes the ability +// to keep track of the time at which the deposit layer was removed. +// Introduced option NOCHANNELDIFFUS to switch off diffusion at channelized +// sites. +// Extended exponent uexp for the U-shape of glacial valleys towards elliptical +// shapes (uexp < 0). +// Fixed problem with combination of options LAKES and DEPOS. +// +// Version 41 2023-02-21 +// Introduced stream-power exponent spexph for hillslopes used in combination +// of the options CHANNEL, SHAREDSP, and SPEXP. +// Introduced exponent uexp for the U-shape of glacial valleys. +// Reorganized default definitions to fix a problem with the combination +// CHANNEL and SHAREDSP. +// +// Version 40 2022-07-06 +// Fixed problem in constructor Grid() in combination with option PRECIP. +// Extended methods for reading so that files written with a different value +// of LAYERS can be read. +// +// Version 39 2022-06-15 +// Changed member spfac to public in order to allow custom stream-power laws. +// Changed sites that contain ice to channelized in order to achive a +// compatibility of the options ICE and CHANNEL. +// Introduced erodibilities kdh and kth for hillslopes used in combination +// of the options CHANNEL and SHAREDSP. +// Removed option RESCALE. +// +// Version 38 2022-04-15 +// Introduced option CHANNEL for the distinction between channel and hillslope +// sites by the flow pattern. +// +// Version 37 2022-01-18 +// Modified implementation of diffusion in such a way that the option +// DIFFSEMIIMPL without DIFFUS and EIGHTNEIGHBORS implements diffusion +// in D8 flow direction. +// Added option DEPOS for switching to fully transport-limited conditons +// in domains of sediment deposition. +// Added option MELTBOUNDARY for melting all ice at the boundaries. +// Fixed a problem in constructing a grid of size 1x1. +// Fixed a problem in reading unrecognized keys. +// Fixed a problem with option MELTRATE. +// +// Version 36 2021-09-08 +// Removed decprecated options SEDIM, CONV and fixed-point iteration for +// method erode(). +// Improved treatment of boundary points for the LFPM. +// +// Version 35 2021-08-19 +// Included LFPM for orographic precipitation. +// +// Version 34 2021-08-03 +// Added option DIFFUSIVITY for spatially variable diffusivity. +// Added option ICEFRAC for using own expressions for the elevation-dependent +// ice production. +// Added option MELTRATE for defining a melting rate instead of the default +// mixing model. +// +// Version 33 2021-03-29 +// Added support for stream-power exponents (SPEXP) in all modes. +// Improved option LAKES for linear decline model. +// Added option LAKEFLOWDIR for improved flow direction in lakes. +// Removed option BOUNDARY. +// Fixed problem in reading array data. +// +// Version 32 2021-03-02 +// Changed interpretation of variable bottom, relative to surface now. +// Added option for writing individual layers of array-valued variables. +// Fixed error in the definition of the keys for erodibilities. +// Added continued reading after finding an unrecognized key in a file. +// +// Version 31 2021-01-19 +// Fixed problem with glacial erodibilities under option NOMELTWATER +// Added file handling for array-valued variables +// Added support for layers with different erodibilities +// +// Version 30 2021-01-08 +// Fixed problems with unrecognized keys. +// Added initializer for variable excav. +// +// Version 29 2020-12-22 +// Fixed problem with option ICE if glacial erosion is not detachment-limited. +// Changed mean of erodibilities for mixed ice/water conditions. +// Shared stream-power model is now default for option ICE. +// Added option NOMELTWATER. +// +// Version 28 2020-12-17 +// Fixed problem with option SHAREDSP. +// Proceeded with work on option ICE. +// +// Version 27 2020-11-27 +// Introduced some technical simplificiations that should not affect the results. +// Proceeded with work on option ICE. +// +// Version 26 2020-10-21 +// Introduced updated methods for obtaining the neighbors. +// Simplified the computation for the transport-limited model and the linear +// decline model. +// Added option LAKES for reducing stream power in lakes. +// Improved information about unrecognized keys when reading files. +// +// Version 25 2020-10-12 +// Fixed problem with lowercase letters in keys. +// +// Version 24 2020-10-06 +// Fixed problem with resizing. +// Introduced option EIGHTNEIGHBORS for a 9-point diffusion stencil. +// Introduced option DIFFSEMIIMPL for semi-implicit diffusion scheme. +// Renamed option SEDIMIMPL TO LINDECL and added option SHAREDSP. +// Replaced option REMOVE by a parameter excav. +// Proceeded with work on option ICE. +// +// Version 23 2020-09-06 +// Added partly implicit scheme for hillslope diffusion for transport-limited +// erosion. Cannot be combined with option REMOVE. +// Removed method diffuse(). +// Renamed variable ac to athr (fluvial threshold). +// Introduced new methods for reading and writing data. +// Reduced writing precision for the discharge with option PRECIP to 4 Bytes. +// Changed meaning of the third argument in method write(). +// Added option ICE (experimental). +// +// Version 22 2020-07-29 +// Removed test part in method erode as this altered the computed sediment +// fluxes +// +// Version 21 2020-06-03 +// Introduced implicit scheme for the transport model of Davy and Lague (2009) +// +// Version 20 2020-04-14 +// Introduced transport-limited erosion TRANSPORT. +// This option is not compatible with the option PRECIP and with a fluvial +// threshold ac > 0. +// Introduced method setErosionFactors(), has to be called explicitly if +// the vales A^theta shall be preset for higher performance. +// Method computeDischarge() is now called automatically by the method erode(). +// Changed definition of class to class template. Use Grid<> instead of Grid. +// +// Version 19 2020-03-09 +// Introduced threshold catchment size ac +// Added new functions for setting parameters +// +// Version 18 2020-02-13 +// Introduced option PRECIP +// +// Version 17 2020-02-03 +// Changed definition of alpha according to Turcotte and Schubert + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace openlem +{ +#define SHAREDSP +#define LAYERS 2 +#ifdef ICE +#define SHAREDSP +#define BUFFERS +#endif + +#ifdef DEPOS +#define SHAREDSP +#define BUFFERS +#endif + +#ifdef SHAREDSP +#define LINDECL +#ifdef CHANNEL +#define BUFFERS +#endif +#endif + +#ifdef LINDECL +#define TRANSPORT +#ifdef LAKES +#define BUFFERS +#endif +#endif + +#ifdef RIGIDITY +#define FLEXURE +#endif + +#ifdef EXTPRECIP +#define PRECIP +#endif + +#define DISTTONEIGH(p,q) sqrt((p.i!=q.i)+(p.j!=q.j)*dx[p.i]*dx[p.i]) + + using namespace std; + + class Keymapentry + { + public: + char keystring[4], varname[16]; + int keyint, offset, size, fsize, grid, readonly, array; + }; + + class Point + { + public: + short int i, j; + + Point ( int i = 0, int j = 0 ) + { + this->i = i; + this->j = j; + } + + void set ( int i, int j ) + { + this->i = i; + this->j = j; + } + + int equals ( int i, int j ) + { + return i == this->i && j == this->j; + } + + int isDiag ( Point p ) + { + return i != p.i && j != p.j; + } + + int isDiag ( int i, int j ) + { + return i != this->i && j != this->j; + } + + int operator== ( Point p ) + { + return i == p.i && j == p.j; + } + + int operator!= ( Point p ) + { + return i != p.i || j != p.j; + } + public: + bool operator== (const Point p) const + { + return i == p.i && j == p.j; + } + + }; + + class Node + { + public: + double h, // surface elevation + l; // water level, >= h + float u; // uplift rate +#ifdef TRANSPORT + double qs; // sediment flux + double qp; // derivative of sediment flux +#endif +#ifdef ICE + double qi; // ice flux + double th; // ice thickness + Point trunk; // donor with the largest catchment size +#endif +#ifdef LAKES +#ifndef ICE + Point trunk; // additional donor in lakes +#endif +#endif +#ifdef ERODIBILITY + float e; // erodibility +#endif +#ifdef DIFFUS + float r; // diffusive erosion rate, does not include fluxes + // in flow direction for transport-limited erosion +#endif +#ifdef DIFFUSIVITY + float diff; // diffusivity +#endif +#ifdef LAYERS + float bottom[LAYERS-1]; // lower boundary of the respective layer + // relative to the surface +#endif +#ifdef FLEXURE + double w; // vertical deflection + float g; // load +#endif +#ifdef RIGIDITY + float alpha; // flexural parameter (4D/(rhom g))^0.25 +#endif +#ifdef PRECIP + float p; // precipitation rate + double q; // discharge +#ifdef EXTPRECIP + float qv; // vapor flux + float qc; // cloud water flux + float ptot; // total precipitation +#endif + +#else + unsigned int q; // discharge (catchment size) +#endif +#ifdef CHANNEL + char channel; // channel marker (1 for channelized sites, 0 else) +#endif + char b, // boundary marker (1 for boundary nodes, 0 else) + m; // marker (only internal use) + Point d; // flow target +#ifdef BUFFERS + double lambdaq; +#endif +#ifdef ICE + double qeff; +#endif + + int drainsTo ( Point p ) + { + return p.i == d.i && p.j == d.j; + } + + int drainsTo ( int i, int j ) + { + return i == d.i && j == d.j; + } + +#ifdef LAYERS + void adjustLayers ( double dh, double t = FLT_MAX ) + { +#ifdef BORROW + if ( dh > 0 || bottom[0] < 0 ) bottom[0] -= dh; + float borrowed = bottom[0] >= 0 ? bottom[0] : 0; + for ( int k = 1; k < LAYERS-1; ++k ) + if ( bottom[k] <= 0 && (bottom[k]-=dh) >= 0 && t > 0 ) bottom[k] = t; +#else + if ( bottom[0] > 0 && dh > 0 ) bottom[0] = 0; + for ( int k = 0; k < LAYERS-1; ++k ) + if ( bottom[k] <= 0 && (bottom[k]-=dh) >= 0 && t > 0 ) bottom[k] = t; +#endif + } + + int getLayer() + { + for ( int i = 0; i < LAYERS-1; ++i ) + if ( bottom[i] < 0 ) return i; + return LAYERS-1; + } +#endif + }; + + class Vect + { + public: + double u1, u2; + + Vect ( double u1 = 0, double u2 = 0 ) + { + this->u1 = u1; + this->u2 = u2; + } + + Vect operator -= ( Vect u ) + { + u1 -= u.u1; + u2 -= u.u2; + return *this; + } + }; + + class Matrix + { + public: + double a11, a12, a21, a22; + + Matrix ( double a11 = 0, double a12 = 0, double a21 = 0, double a22 = 0 ) + { + this->a11 = a11; + this->a12 = a12; + this->a21 = a21; + this->a22 = a22; + } + + Matrix operator -= ( Matrix a ) + { + a11 -= a.a11; + a12 -= a.a12; + a21 -= a.a21; + a22 -= a.a22; + return *this; + } + + Vect operator * ( Vect b ) + { + return Vect ( a11*b.u1+a12*b.u2, a21*b.u1+a22*b.u2 ); + } + + Matrix operator * ( Matrix b ) + { + return Matrix ( a11*b.a11+a12*b.a21, a11*b.a12+a12*b.a22, + a21*b.a11+a22*b.a21, a21*b.a12+a22*b.a22 ); + } + + Matrix inv () + { + double det = a11*a22-a12*a21; + return Matrix ( a22/det, -a12/det, -a21/det, a11/det ); + } + }; + + template + class PointValue + { + public: + Point p; + valuetype d; + + PointValue ( Point p = Point(0,0), valuetype d = 0 ) + { + this->p = p; + this->d = d; + } + }; + + template + class Grid + { + vector > u; // lattice of nodes + Point shift[8]; // relative position of the 8 neighbors + vector imap; // mapping for coarsening + // (only internal use) + vector jmap; // mapping for coarsening + // (only internal use) + unsigned int c; // coarsening level for flexure + // (only internal use) + vector > > w, res, frhs; // for flexure + vector dx; // grid spacing in x-direction for option DX +#ifdef RIGIDITY + vector > > rig; +#else + vector rig; +#endif + double at; // a^theta (only internal use) + + public: + short int m, n; // lattice size m times n + double t, // time + theta, // concavity index + spexp, // exponent in the stream power law + spexph, // exponent in the stream power law at hillslopes + lambda, // sediment deposition parameter in the linear decline model +#ifdef LAYERS + kd[LAYERS], // erodibility for detachment-limited conditions + kt[LAYERS], // erodibility for transport-limited conditions + kdg[LAYERS], // erodibility for detachment-limited conditions + // for glacial erosion + ktg[LAYERS], // erodibility for transport-limited conditions + // for glacial erosion + kdh[LAYERS], // erodibility for detachment-limited conditions + // at hillslopes + kth[LAYERS], // erodibility for transport-limited conditions + // at hillslopes +#else + kd, // erodibility for detachment-limited conditions + kt, // erodibility for transport-limited conditions + kdg, // erodibility for detachment-limited conditions + // for glacial erosion + ktg, // erodibility for transport-limited conditions + // for glacial erosion + kdh, // erodibility for detachment-limited conditions + // at hillslopes + kth, // erodibility for transport-limited conditions + // at hillslopes +#endif + diff, // diffusivity + excav, // factor defining which fraction of the material + // moving from the hillslopes into rivers is + // excavated for detachment-limited erosion +// Option LAKES + lff, // fraction of the actual lake depth attempted to be + // filled with sediments + mld, // minimum lake depth +// Option ICE + ela, // glacial equilibrium line altitude + fsa, // altitude where all precipitation is snow + gwf, gwe, // parameters in the relation + // width = gwf*qi^gwe + ghf, ghe, // not used + ghw, // thickness to width ratio of glaciers + uexp, // exponent for the U-shape of glacial valley profiles + ge, // not used + eps, // weight factor in the geometric mean q^eps*qi^(1-eps) + // used in the stream-power term for qi > q +// Option FLEXURE + alpha, // flexural parameter (4*D/(rhom g))^0.25 + tau, // flexural relaxation time + rhoc, // density ratio crust/mantle + +// Option PRECIP + lc, // length scale of condensation + lf, // length scale of fallout + ll, // length scale of long-range transport + ld, // length scale of dispersion + refheight, // reference elevation + refslope, // reference slope for evaporation + evap, // evaporation fraction at sea level + qin, // influx +// Option CHANNEL + ah; // apparent catchment size for non-channel nodes + unsigned int athr, // threshold catchment size where fluvial erosion starts + a, // shift in catchment size used in the stream-power law + dir; // edge used for influx + vector spfac; // factor in stream power law + map km; + vector diag, upper, lower, right, bottom; + vector rhs; + vector psi; +// Option REDREC + vector seq; + int nrec; + vector prevrow, nextrow, prevcol, nextcol; + class Lake; + vector lakes; // all detected lakes +// double qdev; +// Option LAKES + vector> offset; + +#ifdef DEPOS + int etrue, efalse, dtrue, dfalse, oppos; + int etruesum, efalsesum, dtruesum, dfalsesum, oppossum; +#endif + + Grid ( int m = 1, int n = 1, double theta = 0.5, unsigned int a = 0, + double spexp = 1., double lambda = 1., double diff = 0., + double alpha = 500., double rhoc = 0.85, double tau = 0. ) + { + this->m = this->n = 0; + t = 0.; + setErosionParams(theta,spexp,0,a); + setSharedStreamPowerParams(); + this->lambda = lambda; + this->diff = diff; + excav = 0.; + lff = 1; + mld = 0; +// Default parameters for orographic precipitation +// useful for 100 m grid spacing + lc = lf = ld = 250; + ll = 10000; + refheight = 10; + evap = dir = 0; + qin = 5000; +// Default parameters for glacial erosion + ela = 10; + fsa = 11; + gwf = 1.5; + gwe = 0.3; + ghw = 0.1; + uexp = 0; + eps = 0.25; + + setFlexureParams(alpha,rhoc,tau); + resize ( m, n ); + Point p(0,0); + int k = 0; + for ( p.i = -1; p.i <= 1; ++p.i ) + for ( p.j = -1; p.j <= 1; ++p.j ) + if ( p.i || p.j ) shift[k++] = p; + setKeys(); + } + + void addKey ( const char *key, const char *varname, void *p, unsigned int size, unsigned int fsize, unsigned int grid, int readonly = 0, int array = 1 ) + { + Keymapentry e; + strncpy(e.keystring,key,4); + if ( grid ) + for ( int i = 0; i < 4; ++i ) e.keystring[i] = toupper(e.keystring[i]); + memcpy(&e.keyint,e.keystring,4); + strncpy(e.varname,varname,15); + e.offset = (char *)p - ( grid ? (char *)getNode(0,0) : (char *)&t ); + e.size = size; + e.fsize = fsize; + e.grid = grid; + e.readonly = readonly; + e.array = array; + km[e.keyint] = e; + } + + Keymapentry getKey ( unsigned int keyint ) + { + return km[keyint]; + } + + Keymapentry getKey ( const char *keystring ) + { + unsigned int keyint; + strncpy((char *)&keyint,keystring,4); + return km[keyint]; + } + + void setKeys() + { + addKey ( "t", "t", &t, sizeof(t), 8, 0 ); + addKey ( "c", "theta", &theta, sizeof(theta), 8, 0, 1 ); + addKey ( "thf", "theta", &theta, sizeof(theta), 8, 0 ); + addKey ( "e", "spexp", &spexp, sizeof(spexp), 8, 0, 1 ); + addKey ( "exf", "spexp", &spexp, sizeof(spexp), 8, 0 ); + addKey ( "exh", "spexph", &spexph, sizeof(spexph), 8, 0 ); + addKey ( "l", "lambda", &lambda, sizeof(lambda), 8, 0, 1 ); + addKey ( "sdc", "lambda", &lambda, sizeof(lambda), 8, 0 ); + addKey ( "d", "diff", &diff, sizeof(diff), 8, 0 ); + addKey ( "exc", "excav", &excav, sizeof(excav), 8, 0 ); + addKey ( "lff", "lff", &lff, sizeof(lff), 8, 0 ); + addKey ( "mld", "mld", &mld, sizeof(mld), 8, 0 ); +#ifdef LAYERS + addKey ( "kd", "kd", kd, sizeof(kd[0]), 8, 0, 0, LAYERS ); + addKey ( "kt", "kt", kt, sizeof(kt[0]), 8, 0, 0, LAYERS ); + addKey ( "kdg", "kdg", kdg, sizeof(kdg[0]), 8, 0, 0, LAYERS ); + addKey ( "ktg", "ktg", ktg, sizeof(ktg[0]), 8, 0, 0, LAYERS ); + addKey ( "kdh", "kdh", kdh, sizeof(kdh[0]), 8, 0, 0, LAYERS ); + addKey ( "kth", "kth", kth, sizeof(kth[0]), 8, 0, 0, LAYERS ); +#else + addKey ( "kd", "kd", &kd, sizeof(kd), 8, 0 ); + addKey ( "kt", "kt", &kt, sizeof(kt), 8, 0 ); + addKey ( "kdg", "kdg", &kdg, sizeof(kdg), 8, 0 ); + addKey ( "ktg", "ktg", &ktg, sizeof(ktg), 8, 0 ); + addKey ( "kdh", "kdh", &kdh, sizeof(kdh), 8, 0 ); + addKey ( "kth", "kth", &kth, sizeof(kth), 8, 0 ); +#endif + addKey ( "ela", "ela", &ela, sizeof(ela), 8, 0 ); + addKey ( "fsa", "fsa", &fsa, sizeof(fsa), 8, 0 ); + addKey ( "gwf", "gwf", &gwf, sizeof(gwf), 8, 0 ); + addKey ( "gwe", "gwe", &gwe, sizeof(gwe), 8, 0 ); + addKey ( "ghf", "ghf", &ghf, sizeof(ghf), 8, 0, 1 ); + addKey ( "ghe", "ghe", &ghe, sizeof(ghe), 8, 0, 1 ); + addKey ( "ghw", "ghw", &ghw, sizeof(ghw), 8, 0 ); + addKey ( "uex", "uexp", &uexp, sizeof(uexp), 8, 0 ); + addKey ( "ge", "ge", &ge, sizeof(ge), 8, 0, 1 ); + addKey ( "eps", "eps", &eps, sizeof(eps), 8, 0 ); + addKey ( "p", "alpha", &alpha, sizeof(alpha), 8, 0, 1 ); + addKey ( "alp", "alpha", &alpha, sizeof(alpha), 8, 0 ); + addKey ( "u", "tau", &tau, sizeof(tau), 8, 0, 1 ); + addKey ( "tau", "tau", &tau, sizeof(tau), 8, 0 ); + addKey ( "r", "rhoc", &rhoc, sizeof(rhoc), 8, 0, 1 ); + addKey ( "rho", "rhoc", &rhoc, sizeof(rhoc), 8, 0 ); + addKey ( "lc", "lc", &lc, sizeof(lc), 8, 0 ); + addKey ( "lf", "lf", &lf, sizeof(lf), 8, 0 ); + addKey ( "ll", "ll", &ll, sizeof(ll), 8, 0 ); + addKey ( "ld", "ld", &ld, sizeof(ld), 8, 0 ); + addKey ( "rh", "refheight", &refheight, sizeof(refheight), 8, 0 ); + addKey ( "rs", "refslope", &refslope, sizeof(refslope), 8, 0 ); + addKey ( "ev", "evap", &evap, sizeof(evap), 8, 0 ); + addKey ( "qin", "qin", &qin, sizeof(qin), 8, 0 ); + addKey ( "ah", "ah", &ah, sizeof(ah), 8, 0 ); + addKey ( "dir", "dir", &dir, sizeof(dir), 4, 0 ); + addKey ( "f", "athr", &athr, sizeof(athr), 4, 0, 1 ); + addKey ( "ath", "athr", &athr, sizeof(athr), 4, 0 ); + addKey ( "a", "a", &a, sizeof(a), 4, 0, 1 ); + addKey ( "acr", "a", &a, sizeof(a), 4, 0 ); + + Node *p = getNode(0,0); + addKey ( "H", "h", &p->h, sizeof(p->h), 8, 1 ); + addKey ( "H8", "h", &p->h, sizeof(p->h), 8, 1 ); + addKey ( "H4", "h", &p->h, sizeof(p->h), 4, 1 ); + addKey ( "L", "l", &p->l, sizeof(p->l), 8, 1 ); + addKey ( "L8", "l", &p->l, sizeof(p->l), 8, 1 ); + addKey ( "L4", "l", &p->l, sizeof(p->l), 4, 1 ); + addKey ( "WL", "l", &p->l, sizeof(p->l), 8, 1, 1 ); + addKey ( "WL8", "l", &p->l, sizeof(p->l), 8, 1, 1 ); + addKey ( "WL4", "l", &p->l, sizeof(p->l), 4, 1, 1 ); + addKey ( "U", "u", &p->u, sizeof(p->u), 4, 1 ); +#ifdef TRANSPORT + addKey ( "QS", "qs", &p->qs, sizeof(p->qs), 4, 1 ); +#endif +#ifdef ICE + addKey ( "QI", "qi", &p->qi, sizeof(p->qi), 4, 1 ); + addKey ( "TH", "th", &p->th, sizeof(p->th), 4, 1 ); +#endif +#ifdef ERODIBILITY + addKey ( "E", "e", &p->e, sizeof(p->e), 4, 1 ); +#endif +#ifdef DIFFUS + addKey ( "R", "r", &p->r, sizeof(p->r), 4, 1 ); +#endif +#ifdef DIFFUSIVITY + addKey ( "DIF", "diff", &p->diff, sizeof(p->diff), 4, 1 ); +#endif +#ifdef LAYERS + addKey ( "BOT", "bottom", p->bottom, sizeof(p->bottom[0]), 4, 1, 0, LAYERS-1 ); +#endif +#ifdef FLEXURE + addKey ( "W", "w", &p->w, sizeof(p->w), 8, 1 ); + addKey ( "W8", "w", &p->w, sizeof(p->w), 8, 1 ); + addKey ( "W4", "w", &p->w, sizeof(p->w), 4, 1 ); + addKey ( "G", "g", &p->g, sizeof(p->g), 4, 1 ); + addKey ( "LOA", "g", &p->g, sizeof(p->g), 4, 1 ); +#endif +#ifdef RIGIDITY + addKey ( "X", "alpha", &p->alpha, sizeof(p->alpha), 4, 1, 1 ); + addKey ( "ALP", "alpha", &p->alpha, sizeof(p->alpha), 4, 1 ); +#endif +#ifdef PRECIP + addKey ( "P", "p", &p->p, sizeof(p->p), 4, 1 ); +#endif +#ifdef EXTPRECIP + addKey ( "QV", "qv", &p->qv, sizeof(p->qv), 4, 1 ); + addKey ( "QC", "qc", &p->qc, sizeof(p->qc), 4, 1 ); + addKey ( "PT", "ptot", &p->ptot, sizeof(p->ptot), 4, 1 ); +#endif + addKey ( "Q", "q", &p->q, sizeof(p->q), 4, 1 ); +#ifdef CHANNEL + addKey ( "CHA", "channel", &p->channel, sizeof(p->channel), 1, 1 ); +#endif + addKey ( "B", "b", &p->b, sizeof(p->b), 1, 1 ); + addKey ( "D", "d", &p->d, sizeof(p->d), sizeof(p->d), 1 ); + addKey ( "S", "s", &p, 0, sizeof(float), 1 ); + } + + void setSharedStreamPowerParams ( int layer = 0, double kd = 2, double kt = 2, + double kdg = 10, double ktg = 1e100, + double kdh = 6, double kth = 6 ) + { +#ifdef LAYERS + this->kd[layer] = kd; + this->kt[layer] = kt; + this->kdg[layer] = kdg; + this->ktg[layer] = ktg; + this->kdh[layer] = kdh; + this->kth[layer] = kth; +#else + this->kd = kd; + this->kt = kt; + this->kdg = kdg; + this->ktg = ktg; + this->kdh = kdh; + this->kth = kth; +#endif + } + + void setErosionParams ( double theta = 0.5, double spexp = 1., + unsigned int athr = 0, unsigned int a = 0 ) + { + this->theta = theta; + this->spexp = spexp; + this->spexph = 1; + this->athr = athr; + this->a = a; + this->ah = 9; + } + + void setFlexureParams ( double alpha = 500., double rhoc = 0.85, double tau = 0. ) + { + this->alpha = alpha; + this->rhoc = rhoc; + this->tau = tau; + } + + void setPrecipParams ( double lc, double lf, double ll, + double ld, double refheight, double refslope, + double evap, double qin, int dir = 0 ) + { + this->lc = lc; + this->lf = lf; + this->ll = ll; + this->ld = ld; + this->refheight = refheight; + this->refslope = refslope; + this->evap = evap; + setInflux ( qin, dir ); + } + + void setInflux ( double qin, int dir = 0 ) + { + this->qin = qin; + this->dir = dir; + } + + void resize ( int m, int n ) + { + if ( this->m == m && this->n == n ) return; + this->m = m; + this->n = n; + u.resize(m); + prevrow.resize(m); + prevrow[0] = m-1; + for ( int i = 1; i < m; ++i ) prevrow[i] = i-1; + nextrow.resize(m); + for ( int i = 0; i < m-1; ++i ) nextrow[i] = i+1; + nextrow[m-1] = 0; + prevcol.resize(n); + prevcol[0] = n-1; + for ( int j = 1; j < n; ++j ) prevcol[j] = j-1; + nextcol.resize(n); + for ( int j = 0; j < n-1; ++j ) nextcol[j] = j+1; + nextcol[n-1] = 0; +#ifdef DX + dx.resize(m,1.); +#endif +#ifdef REDREC + seq.clear(); + seq.reserve(m*n); +#endif + for ( int i = 0; i < m; ++i ) + { + u[i].resize(n); + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + pn->h = pn->l = pn->u = pn->b = 0; + pn->d.set(i,j); +#ifdef TRANSPORT + pn->qs = pn->qp = 0; +#endif +#ifdef CHANNEL + pn->channel = 0; +#endif +#ifdef PRECIP + pn->p = 1; +#endif +#ifdef ERODIBILITY + pn->e = 1.; +#endif +#ifdef DIFFUSIVITY + pn->diff = diff; +#endif +#ifdef RIGIDITY + pn->alpha = alpha; +#endif +#ifdef ICE + pn->qi = pn->th = 0; +#endif +#ifdef FLEXURE + pn->w = pn->g = 0; +#endif +#ifdef REDREC + seq.push_back(Point(i,j)); +#endif + } + } +#ifdef FLEXURE + w.clear(); + res.clear(); + frhs.clear(); + rig.clear(); + do + { + w.push_back(*new vector >); + w.back().resize(m); + for ( int i = 0; i < m; ++i ) w.back()[i].resize(n); + res.push_back(*new vector >); + res.back().resize(m); + for ( int i = 0; i < m; ++i ) res.back()[i].resize(n); + frhs.push_back(*new vector >); + frhs.back().resize(m); + for ( int i = 0; i < m; ++i ) frhs.back()[i].resize(n); +#ifdef RIGIDITY + rig.push_back(*new vector >); + rig.back().resize(m); + for ( int i = 0; i < m; ++i ) rig.back()[i].resize(n); +#else // NOT RIGIDITY + rig.push_back(0.); +#endif // NOT RIGIDITY + } + while ( m%2 == 0 && n%2 == 0 && (m/=2)>=2 && (n/=2)>=2 ); +#endif // FLEXURE +#ifdef LAKES + offset.resize(m*n); + int k = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + offset[k].p = Point(i,j); + int di = i <= m/2 ? i : m-i; + int dj = j <= n/2 ? j : n-j; + offset[k++].d = di*di+dj*dj; + } + offset[0].d = -1; + sort(offset); +#endif +#ifdef PRECIP + n = this->n >= 2 ? this->n : 2; + if ( this->m > n ) n = this->m; + diag.resize(n); + upper.resize(n-1); + lower.resize(n-1); + right.resize(n-2); + bottom.resize(n-2); + rhs.resize(n); + rhs.resize(n); + psi.resize(n); +#endif // PRECIP + + } + +#ifndef PRECIP + void setErosionFactors() + { + spfac.resize(m*n); + at = pow(a,theta); + for ( int q = 0; q < spfac.size(); ++q ) + spfac[q] = q >= athr ? pow(q,theta)+at : 0.; + } +#endif + + Node *getNode ( int i, int j ) + { + return &u[i][j]; + } + + Node *getNode ( Point p ) + { + return &u[p.i][p.j]; + } + + Node *getNode ( Point *p ) + { + return &u[p->i][p->j]; + } + + Node *getNode ( vector::iterator p ) + { + return &u[p->i][p->j]; + } + + Node *getNodeP ( int i, int j ) + { + return &u[(i+m)%m][(j+n)%n]; + } + + vector &operator [] ( int i ) + { + return u[i]; + } + + Node &operator [] ( Point p ) + { + return u[p.i][p.j]; + } + + Node &operator [] ( Point *p ) + { + return u[p->i][p->j]; + } + + Node &operator [] ( vector::iterator p ) + { + return u[p->i][p->j]; + } + + Point getNeighbor ( Point p, int i ) + { + return Point ( (p.i+shift[i].i+m)%m, (p.j+shift[i].j+n)%n ); + } + + vector getNeighbors ( Point p ) +// Returns a vector of the 8 nearest and second nearest neighbors + { + vector neigh = { Point(prevrow[p.i],prevcol[p.j]), + Point(prevrow[p.i],p.j), + Point(prevrow[p.i],nextcol[p.j]), + Point(p.i,prevcol[p.j]), + Point(p.i,nextcol[p.j]), + Point(nextrow[p.i],prevcol[p.j]), + Point(nextrow[p.i],p.j), + Point(nextrow[p.i],nextcol[p.j]) + }; + return neigh; + } + + vector getNearestNeighbors ( Point p ) +// Returns a vector of the 4 nearest neighbors + { + vector neigh = { Point(prevrow[p.i],p.j), + Point(p.i,prevcol[p.j]), + Point(p.i,nextcol[p.j]), + Point(nextrow[p.i],p.j) + }; + return neigh; + } + + int computeDistSquare ( Point p, Point q ) + { + int di = abs(p.i-q.i); + if ( 2*di > m ) di = m-di; + int dj = abs(p.j-q.j); + if ( 2*dj > n ) dj = n-dj; + return di*di+dj*dj; + } + +#ifdef REDREC +#ifdef DOWN + void setSeq() +// Defines a sequence of points in downstream order. Boundary points are +// not included and marked with a positive value. + { + int nb = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) nb += getNode(i,j)->m = getNode(i,j)->b; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) ++getNode(getNode(i,j)->d)->m; + seq.clear(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Point p(i,j); + Node *pn = getNode(p); + while ( pn->m == 0 ) + { + seq.push_back(p); + --pn->m; + /* + if ( pn->d == p ) + { + printf ( "Fatal error\n" ); + exit(0); + } + */ + p = pn->d; + --(pn=getNode(p))->m; + } + } + for ( vector::iterator p = seq.begin(); p != seq.end(); ++p ) + getNode(p)->m = 0; + if ( seq.size()+nb != m*n ) + { + printf ( "Fatal error: Wrong size of sequence %i\n", seq.size() ); + exit ( 0 ); + } + /* + printf ( "%i\n", seq.size() ); + vector num(20,0); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) ++num[getNode(i,j)->m+10]; + for ( int i = 0; i < num.size(); ++i ) + printf ( "%i %i\n", i-10, num[i] ); + */ + } + +#else + + void setSeq() + { + seq.clear(); + vector::iterator s = seq.begin(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->b ) + { + seq.push_back(Point(i,j)); + for ( ; s != seq.end(); ++s ) + { + vector neigh = getNeighbors(*s); + for ( vector::iterator p = neigh.begin(); p != neigh.end(); ++p ) + if ( getNode(p)->drainsTo(*s) ) + seq.push_back(*p); +#ifdef CHANNEL + else + getNode(s)->channel -= ( getNode(s)->h >= getNode(p)->h ); +#endif +#ifdef LAKES +#ifndef LAKEFLOWDIR + if ( getNode(*s)->trunk != *s ) seq.push_back(getNode(*s)->trunk); +#endif +#endif + } + } + if ( seq.size() != m*n ) + { + printf ( "Fatal error: Wrong size of sequence %i\n", seq.size() ); + exit ( 0 ); + } + } + +#endif +#endif + + private: + int computeWaterLevel ( Point p ) + { + int nc = 0; + double l = 0.; + if ( getNode(p)->m ) return nc; + getNode(p)->m = 1; + l = getNode(p)->h; + Point dest = getNode(p)->d; + if ( dest != p ) nc += computeWaterLevel(dest); + if ( getNode(dest)->l > l ) l = getNode(dest)->l; + if ( l != getNode(p)->l ) ++nc; + getNode(p)->l = l; + return nc; + } + + public: + int computeWaterLevel ( ) +// Recomputes the water level using given flow directions and returns the number +// of sites where the water level has changed. Calling this function makes sense +// after changing surface elevations by any other process than fluvial erosion +// before recomputing the flow directions. If the surface has been changed only +// by fluvial erosion, calling this method has no effect. + { + int nc = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + nc += computeWaterLevel(Point(i,j)); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + return nc; + } + +#ifdef LAKEFLOWDIR + unsigned int computeFlowDirection ( Point p ) + { +// Computes the flow direction of a node. Leaves it unchanged if no downward +// direction was found or if the steepest downward direction results in an +// invalid flow pattern. + Node *pn = getNode(p); + if ( pn->b ) return 0; + Point pd = pn->d; + double maxs = 0.; + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + { + double s = pn->h - getNode(u)->h; + if ( p.isDiag(*u) ) s *= sqrt(0.5); + if ( s > maxs ) + { + pd = *u; + maxs = s; + } + } + if ( pd == pn->d ) return 0; + Point d = pd; + Node *dn; + while ( !(dn=getNode(d))->b ) + { + if ( d == p ) + { +// pn->q = 1; +// printf ( "Loop %i %i\n", p.i, p.j ); + return 0; + } + Node *dn = getNode(d); + if ( dn->l > pn->l ) + { +// pn->q = 2; +// printf ( "Increasing water level %i %i %g %i %i %g\n", p.i, p.j, pn->l, d.i, d.j, dn->l ); + return 0; + } + if ( dn->l < pn->l ) + { +// if ( pn->l > pn->h ) +// printf ( "Success %i %i %g %i %i %g\n", p.i, p.j, pn->l, d.i, d.j, dn->l ); + pn->d = pd; + return 1; + } + d = dn->d; + } +// if ( pd != d ) +// printf ( "Boundary %i %i %g %i %i %g\n", p.i, p.j, pn->l, d.i, d.j, dn->l ); + pn->d = pd; + return 1; + } + +#else // NOT LAKEFLOWDIR + unsigned int computeFlowDirection ( Point p ) + { +// Computes the flow direction of a node. Leaves it unchanged if no downward +// direction was found. + Node *pn = getNode(p); + if ( pn->b ) return 0; + Point pd = pn->d; + double maxs = 0.; + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + { + double s = pn->l - getNode(u)->l; + if ( p.isDiag(*u) ) s *= sqrt(0.5); + if ( s > maxs ) + { + pn->d = *u; + maxs = s; + } + } + return pd != pn->d; + } +#endif // NOT LAKEFLOWDIR + + unsigned int computeFlowDirection() +// Computes the flow direction of all nodes based on the water level (not +// surface elevation). Leaves those nodes unchanged where no downward +// direction was found. Returns the number of nodes where the flow direction +// has changed. + { + int nchanges = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + nchanges += computeFlowDirection(Point(i,j)); + return nchanges; + } + +// private: +#ifdef ICE + void propagateUpstream ( Point p, Point pb ) + { + Point pr = pb; + for ( int i = 0; ; ++i ) + { + int di = 2*abs(p.i-pr.i); + if ( di > m ) di = 2*m-di; + int dj = 2*abs(p.j-pr.j); + if ( dj > n ) dj = 2*n-dj; + double dist; + if ( ( dist = (di*di+dj*dj)/(getNode(pr)->th*getNode(pr)->th) ) <= 1 ) + { + if ( getNode(pb)->qi > getNode(p)->qi ) + { + getNode(p)->qi = getNode(pb)->qi; + getNode(p)->qp = dist; + } + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(*u)->drainsTo(p) ) + propagateUpstream(*u,pb); + return; + } + Point pnew = getNode(pr)->trunk; + if ( pnew == pr ) return; + pr = pnew; + } + } + + double computeIceFlux ( Point p ) + { + Node *pn = getNode(p); + if ( pn->m ) return pn->qi; + pn->m = 1; +#ifdef ICEFRAC + double f = (ICEFRAC); +#else + double f = (pn->h+pn->th-ela)/(fsa-ela); +#endif +#ifdef MELTRATE + if ( f < 0 ) f = 0; +#endif +#ifdef PRECIP + pn->q = pn->p; + pn->qi = f*pn->p; + double qin, qmax = 0; +#else + ++pn->q; + pn->qi += f; + int qin, qmax = 0; +#endif +#ifdef MELTBOUNDARY + if ( pn->b ) return pn->qi = pn->th = 0; +#endif + vector neigh = getNeighbors(p); + pn->trunk = p; + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(u)->drainsTo(p) ) + { + pn->qi += computeIceFlux(*u); + pn->q += qin = getNode(u)->q; + if ( qin > qmax ) + { + qmax = qin; + pn->trunk = *u; + } + } + if ( pn->qi > 0. ) + { + if ( pn->qi > pn->q ) pn->qi = pn->q; +#ifdef MELTRATE + f = (MELTRATE); +//printf ( "th %i %i %e\n", p.i, p.j, pn->qi ); + if ( f > 0 && ( pn->qi -= f*gwf*pow(pn->qi,gwe) ) < 0 ) pn->qi = 0; +#endif +// Set thickness to glacier width for the moment + pn->th = gwf*pow(pn->qi,gwe); +//printf ( "th %i %i %e\n", p.i, p.j, pn->qi ); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(u)->drainsTo(p) && *u != pn->trunk ) +// if ( getNode(u)->drainsTo(p) && getNode(u)->q < qmax ) + propagateUpstream ( *u, p ); + } +#ifdef MELTRATE + else pn->qi = pn->th = 0; +#else + else pn->th = 0; +#endif + return pn->qi; + } + +// public: + void computeIceFlux() + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->qp = getNode(i,j)->qi = getNode(i,j)->q = getNode(i,j)->m = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + computeIceFlux ( Point(i,j) ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + if ( getNode(i,j)->qi < 0 ) getNode(i,j)->qi = 0; + getNode(i,j)->m = 0; + } + } + + void computeIceThickness ( Point p ) + { + Node *pn = getNode(p); + if ( pn->m ) return; + pn->m = 1; + pn->th *= ghw; + Point dest = pn->d; + Node *pd = getNode(dest); + if ( dest != p ) + { + computeIceThickness(dest); +// not sure whether this is completely safe + if ( pd->qi == pn->qi && pd->th > pn->th ) pn->th = pd->th; + } + } + + double addIceLayer ( Point p ) + { + double v = 0., dh; + Node *pn = getNode(p); + if ( pn->m ) return 0; + pn->m = 1; + if ( uexp > 0 && pn->qp > 0 ) pn->th *= 1-pow(pn->qp,uexp/2); + if ( uexp < 0 && pn->qp > 0 ) pn->th *= sqrt(1-pn->qp); + Point dest = pn->d; + Node *pd = getNode(dest); + if ( dest != p ) + { + v += addIceLayer(dest); + if ( pd->th > 0 && ( dh = pd->h-pn->h-pn->th ) > 0 ) pn->th += dh; + } + v += pn->th; + pn->h += pn->th; + return v; + } + + double addIceLayer ( ) + { + double v = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + computeIceThickness ( Point(i,j) ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + v += addIceLayer ( Point(i,j) ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + return v; + } + + void removeIceLayer ( Point p ) + { + Node *pn = getNode(p); + if ( pn->m ) return; + pn->m = 1; + pn->h -= pn->th; + pn->l = pn->h; + Point dest = pn->d; + if ( dest != p ) + { + removeIceLayer(dest); + if ( getNode(dest)->l > pn->l ) pn->l = getNode(dest)->l; + } + } + + void removeIceLayer ( ) + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + removeIceLayer ( Point(i,j) ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + } +#endif // ICE + +#ifdef PRECIP + double computeDischarge ( Point p ) +#else + int computeDischarge ( Point p ) +#endif + { +// Computes the discharge of a node + Node *pn = getNode(p); + if ( pn->m ) return pn->q; + pn->m = 1; +#ifdef PRECIP + pn->q += pn->p; +#else + ++(pn->q); +#endif + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(u)->drainsTo(p) ) + pn->q += computeDischarge(*u); + return pn->q; + } + +// public: + void computeDischarge() + { +// Computes all discharges +#ifdef REDREC + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) +#ifdef PRECIP + getNode(i,j)->q += getNode(i,j)->p; +#else + getNode(i,j)->q = 1; +#endif + setSeq(); +#ifdef DOWN + for ( vector::iterator p = seq.begin(); p != seq.end(); ++p ) + getNode(getNode(*p)->d)->q += getNode(*p)->q; +#else + for ( vector::reverse_iterator p = seq.rbegin(); p != seq.rend(); ++p ) + if ( !getNode(*p)->b ) + getNode(getNode(*p)->d)->q += getNode(*p)->q; +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; +#ifndef PRECIP + int q = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->b ) q += getNode(i,j)->q; + if ( q != m*n ) + { + printf ( "Error: Total discharge is %i, should be %i\n", q, m*n ); + exit ( -1 ); + } +#endif +#else + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + pn->q = pn->m = 0; + } + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + computeDischarge ( Point(i,j) ); +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + } + +#ifdef PRECIP + double computeFluxes ( Point p, double dt = 0. ) +#else + int computeFluxes ( Point p, double dt = 0. ) +#endif + { +// Computes the discharge of a node + Node *pn = getNode(p); + if ( pn->m ) return pn->q; + pn->m = 1; + Point dest = pn->d; + Node *destn = getNode(dest); +#ifdef TRANSPORT + double qssum = 0., qpsum = 0.; +#endif + vector neigh = getNeighbors(p); +#ifdef LAKES +#ifndef LAKEFLOWDIR + if ( pn->trunk != p ) neigh.push_back(pn->trunk); +#endif +#endif + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + { + if ( getNode(u)->drainsTo(p) ) + { + if ( !getNode(u)->m ) ++nrec; +#ifdef REDREC +#ifdef DOWN + pn->q += getNode(u)->q; +#endif +#else +// pn->q += computeFluxes(*u,dt); + computeFluxes(*u,dt); +#endif +#ifdef TRANSPORT + qssum += getNode(u)->qs; + qpsum += getNode(u)->qp; +#endif +#ifdef CHANNEL +#if !defined(REDREC) || defined(DOWN) +// pn->channel += 10*getNode(u)->channel; +#endif +#endif + } +#ifdef CHANNEL +#if !defined(REDREC) || defined(DOWN) + else + pn->channel -= ( pn->h >= getNode(u)->h ); +#endif +#endif + } +#ifdef CHANNEL + pn->channel = pn->channel > 0 && !pn->b; +#endif +#ifdef TRANSPORT +// if ( !pn->b ) +// { +// if ( fabs((qssum-pn->qs)/qssum) > qdev ) qdev = fabs((qssum-pn->qs)/qssum); +// if ( fabs((qpsum-pn->qp)/qpsum) > qdev ) qdev = fabs((qpsum-pn->qp)/qpsum); +// } + qssum = pn->qs; + qpsum = pn->qp; + if ( !pn->b && dt > 0. ) + { + double alpha = 1-dt*qpsum; + double q = pn->q; + double spf = 1e-10; + double u = pn->u; +#ifdef DIFFUSIVITY + double diff = pn->diff; +#else + double diff = this->diff; +#endif +#ifdef DIFFUS + u += pn->r; +#ifdef EIGHTNEIGHBORS + diff *= sqrt(2)-1; +#endif // EIGHTNEIGHBORS +#endif // DIFFUS + +#ifdef CHANNEL +#ifdef NOCHANNELDIFFUS + if ( pn->channel ) diff = 0; +#endif +#endif + +#ifdef ICE + double ice = pn->qi/pn->q; + if ( ice > 1 ) + { + q = pn->qi*pow(ice,-eps); + ice = 1; + } +#ifdef NOMELTWATER + else if ( ice > 0 ) + { + q = pn->qi; + ice = 1; + } +#endif // NOMELTWATER +// else if ( ice < 0 ) ice = 0; + pn->qeff = q; + +#endif // ICE + +#ifdef LAKES +#ifndef LINDECL + if ( pn->l <= pn->h || dt*(u+qssum) >= lff*(pn->l-pn->h)*(1-dt*qpsum) ) +#endif // NOT LINDECL +#endif // LAKES + +#ifdef CHANNEL +#ifdef ICE + if ( !pn->channel && ice == 0 ) spf = 1; + else +#else +#ifdef SHAREDSP + if ( !pn->channel ) spf = 1; + else +#else if ( !pn->channel ) spf = pow(ah,theta); + else +#endif +#endif +#endif +#ifdef ICE + spf = q >= athr ? pow(q,theta)+at : 0.; +#else // NOT ICE +#ifdef PRECIP + spf = q >= athr ? pow(q,theta)+at : 0.; +#else // NOT PRECIP + spf = spfac.size() ? spfac[pn->q] : ( q >= athr ? pow(q,theta)+at : 0. ); +#endif // NOT PRECIP +#endif // NOT ICE + +#ifdef LINDECL + double lambda = this->lambda; +#ifdef SHAREDSP + +#ifdef LAYERS + int layer = pn->getLayer(); + double kd = this->kd[layer]; + double kt = this->kt[layer]; + double kdg = this->kdg[layer]; + double ktg = this->ktg[layer]; + double kdh = this->kdh[layer]; + double kth = this->kth[layer]; +#else // NOT LAYERS + double kd = this->kd; + double kt = this->kt; +#endif // NOT LAYERS + +#ifdef CHANNEL +#ifdef ICE + if ( !pn->channel && ice == 0 ) + { + kd = kdh; + kt = kth; + } +#else + if ( !pn->channel ) + { + kd = kdh; + kt = kth; + double spexp = pn->channel ? spexph : this->spexp; + } +#endif +#endif // CHANNEL + + lambda = kd/kt; + +#ifdef ICE + if ( ice > 0 ) + if ( ice >= 1 ) + { + kd = kdg; + lambda = kdg/ktg; + } + else + { + double tmp1 = (1-ice)*kt/kd; + double tmp2 = ice*ktg/kdg; + lambda = 1 / ( tmp1 > tmp2 ? tmp1 : tmp2 ); + kd = pow(1-ice,theta*spexp)*kd + pow(ice,theta*spexp)*kdg; + } +#endif // ICE +#endif // SHAREDSP +#endif // LINDECL + + double f = spf; + double df = p.isDiag(dest) ? sqrt(0.5) : 1.; +#ifdef ERODIBILITY + f *= pn->e; +#endif // ERODIBILITY + +#ifdef SPEXP + f *= pow(spf*fabs(pn->h-destn->h)*df,spexp-1.); +#endif // SPEXP + +#ifdef LINDECL +#ifdef LAKES + if ( pn->l > pn->h ) + { + double me = lff*(pn->l-pn->h)/dt-u; + if ( me > 0 ) + { + double qstot = qssum+lff*(pn->l-pn->h)*qpsum-me; +// double lambdanew = qstot > 0 ? q*me/qstot : 1e100; +// if ( lambdanew > lambda ) lambda = lambdanew; + lambda = qstot > 0 ? q*me/qstot : 1e100; + f = 0; + } + } +#endif // LAKES + +#ifdef SHAREDSP + f *= kd; +#else // NOT SHAREDSP + f *= 1+lambda; +#endif // NOT SHAREDSP + +#ifdef DIFFSEMIIMPL +#ifdef EIGHTNEIGHBORS + f += lambda*diff/q; +#else + f += lambda*diff/q*df; +#endif +#endif // DIFFSEMIIMPL + + f *= df; + double den = alpha*lambda/q+1+dt*f; + pn->qs = (alpha*(f*(pn->h-destn->h)-u)+(u+qssum)*(1+dt*f))/den; + +#ifdef DEPOS +#ifdef LAKES + if ( pn->l <= pn->h ) +#endif + if ( lambda/q*pn->qs > f*(pn->h+dt*u-destn->h) ) + { + f *= 1e100; + lambda *= 1e100; + den = alpha*lambda/q+1+dt*f; + pn->qs = (alpha*(f*(pn->h-destn->h)-u)+(u+qssum)*(1+dt*f))/den; + } +#endif + pn->qp = -alpha*f/den; +#ifdef BUFFERS + pn->lambdaq = lambda/q; +#endif // BUFFERS + +#else // NO LINDECL + + f *= q; + +#ifdef DIFFSEMIIMPL +#ifdef EIGHTNEIGHBORS + f += diff; +#else + f += diff*df; +#endif +#endif // DIFFSEMIIMPL + + f *= df; + double den = dt+alpha/f; + pn->qp = -alpha/den; + pn->qs = (dt*(u+qssum)+alpha*(pn->h-destn->h))/den; +#endif // NO LINDECL + pn->l = f; +#ifndef DOWN + destn->qs += pn->qs; + if ( !destn->b ) destn->qp += pn->qp; +#endif + } +#endif // TRANSPORT + if ( !pn->b ) + { + destn->q += pn->q; +#ifdef CHANNEL + destn->channel += 10*pn->channel; +#endif + } + return pn->q; + } + + int computeFluxes ( int i, int j, double dt = 0. ) + { + Point p; + p.set(i,j); + return computeFluxes ( p, dt ); + } + +// public: + void computeFluxes ( double dt = 0. ) + { +// Computes all discharges and checks for a consistent drainage pattern +#ifdef ICE + computeIceFlux(); + addIceLayer(); +#endif +#ifdef LAKES +#ifndef LAKEFLOWDIR + findLakes(); +#endif +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + pn->q = pn->m = 0; +#ifdef PRECIP + pn->q = pn->p; +#else + pn->q = 1; +#endif +#ifdef TRANSPORT + pn->qs = pn->qp = 0.; +#endif +#ifdef CHANNEL + pn->channel = pn->l == pn->h ? 2 : 10; +#endif +#ifdef DIFFUS + pn->r = 0.; +#endif + } +#ifdef DIFFUS + double diff = this->diff; +#ifdef EIGHTNEIGHBORS + diff *= sqrt(2)-1; +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Point p(i,j); + Node *pn = getNode(p); +#ifdef DIFFUSIVITY + diff = pn->diff; +#ifdef EIGHTNEIGHBORS + diff *= sqrt(2)-1; +#endif // EIGHTNEIGHBORS +#endif // DIFFUSIVITY + +#ifdef CHANNEL +#ifdef NOCHANNELDIFFUS + if ( pn->channel ) diff = 0; +#endif +#endif + +#ifdef EIGHTNEIGHBORS + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) +#ifdef TRANSPORT +#ifdef DIFFSEMIIMPL + if ( !pn->drainsTo(*u) ) +#endif +#endif + { + double dh = pn->h-getNode(u)->h; + if ( dh > 0. ) + { + double flux = diff*dh; + if ( u->isDiag(p) ) flux *= sqrt(0.5); + pn->r -= flux; + getNode(u)->r += flux; + } + } +#else // NOT EIGHTNEIGHBORS + for ( int k = i-1; k <= i+1; k += 2 ) + { + int kper = (k+m)%m; + double dh = pn->h-getNode(kper,j)->h; + if ( dh > 0. ) +#ifdef TRANSPORT +#ifdef DIFFSEMIIMPL + if ( kper != pn->d.i ) +#endif +#endif + { + double flux = diff*dh; + pn->r -= flux; + getNode(kper,j)->r += flux; + } + } + for ( int k = j-1; k <= j+1; k += 2 ) + { + int kper = (k+n)%n; + Node *q = getNode(i,kper); + double dh = pn->h-getNode(i,kper)->h; + if ( dh > 0. ) +#ifdef TRANSPORT +#ifdef DIFFSEMIIMPL + if ( kper != pn->d.j ) +#endif +#endif + { + double flux = diff*dh; + pn->r -= flux; + getNode(i,kper)->r += flux; + } + } +#endif // NOT EIGHTNEIGHBORS + } +#endif // DIFFUS + + nrec = 0; +// qdev = 0; +#ifdef REDREC + setSeq(); + +#ifdef DOWN + for ( vector::iterator p = seq.begin(); p != seq.end(); ++p ) +#else + for ( vector::reverse_iterator p = seq.rbegin(); p != seq.rend(); ++p ) +#endif + { + computeFluxes ( *p, dt ); + /* + Node *pn = getNode(*p); + if ( pn->m ) + { + printf ( "Fatal error\n" ); + exit(0); + } + int tmp = computeFluxes ( *p, dt ); + if ( pn->b ) q += tmp; + */ + } +// printf ( "qdev = %e\n", qdev ); +#else + int q = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + int tmp = computeFluxes ( i, j, dt ); + if ( getNode(i,j)->b ) q += tmp; + } +#ifndef PRECIP + if ( q != m*n ) + { + printf ( "Error: Total discharge is %i, should be %i\n", q, m*n ); + exit ( -1 ); + } +#endif +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; + } + + private: + double erode ( double dt, Point p ) +// Performs a timestep of length dt for a given point and returns +// the absolute change in elevation + { + double e, maxc = 0.; + Node *pn = getNode(p); + Point dest = pn->d; +#ifndef REDREC + if ( pn->m ) return 0; + pn->m = 1; + if ( pn->b ) return 0; + if ( dest == p ) + { + printf ( "Error: Detected node without valid flow target: %i %i\n", + p.i, p.j ); + exit ( -1 ); + } +#ifndef PRECIP + if ( pn->q == 0 ) + { + printf ( "Error: Detected node without discharge: %i %i\n", + p.i, p.j ); + exit ( -1 ); + } +#endif + if ( ( e = erode(dt,dest) ) > maxc ) maxc = e; +#endif + double hold = pn->h; + double q = pn->q; + double spf; + double u = pn->u; + +#ifdef TRANSPORT + +#ifdef DIFFUS + u += pn->r; +#endif + + pn->qs += pn->qp*getNode(dest)->qp; + +#ifdef LINDECL + +#ifdef SHAREDSP + +#ifdef BUFFERS + double lambdaq = pn->lambdaq; +#ifdef ICE + q = pn->qeff; +#endif // ICE + +#else // NO BUFFERS + +#ifdef LAYERS + int layer = pn->getLayer(); + double kd = this->kd[layer]; + double kt = this->kt[layer]; + double kdg = this->kdg[layer]; + double ktg = this->ktg[layer]; +#else // NO LAYERS + double kd = this->kd; +#endif // NO LAYERS + double lambda = kd/kt; + double lambdaq = lambda/q; +#endif // NO BUFFERS + +#else // NO SHAREDSP + +#ifdef BUFFERS + double lambdaq = pn->lambdaq; +#else // NO BUFFERS + double lambdaq = lambda/q; +#endif + +#endif // NO SHAREDSP + double hnew = (pn->h+dt*(pn->l*getNode(dest)->h+lambdaq*pn->qs+u))/(1+dt*pn->l); + +#ifdef DEPOS + if ( pn->qs < 0 ) ++oppos; + if ( lambdaq > 1e10 ) + if ( hnew > pn->h+u*dt ) ++dtrue; + else + { + ++efalse; +// printf ( "%i %e\n", pn->q, (pn->h-hnew)/dt+u ); + } + else if ( hnew > pn->h+u*dt ) ++dfalse; + else ++etrue; +#endif + +#else // NO LINDECL + double hnew = getNode(dest)->h+pn->qs/pn->l; +#endif // NO LINDECL + + pn->qp = hnew-pn->h; + pn->h = hnew; + pn->l = pn->h >= getNode(dest)->l ? pn->h : getNode(dest)->l; + +#else // Detachment-limited + +#ifdef PRECIP + spf = q >= athr ? pow(q,theta)+at : 0.; +#else + spf = spfac.size() ? spfac[pn->q] : ( q >= athr ? pow(q,theta)+at : 0. ); +#endif + double f = spf; + double df = p.isDiag(dest) ? sqrt(0.5) : 1.; +#ifdef ERODIBILITY + f *= pn->e; +#endif + +#ifdef SPEXP + if ( pn->l <= getNode(dest)->l ) f = 0; + else f *= pow(spf*(pn->l-getNode(dest)->l)*df,spexp-1.); +#endif + +// e = 0.; + +#ifdef DIFFUS + u += f > 0 && pn->r > 0 ? (1-excav)*pn->r : pn->r; +#endif + + f *= df; + pn->l = (pn->h+dt*(u+f*getNode(dest)->l))/(1.+dt*f); + if ( pn->l < getNode(dest)->l ) + { + pn->l = getNode(dest)->l; + pn->h += dt*u; + } + else + pn->h = pn->l; +#endif + double dh = pn->h-hold; + if ( ( e = fabs(dh) ) > maxc ) maxc = e; + return maxc; + } + + public: + double erode ( double dt ) +// Performs a timestep of length dt for the entire grid and returns the +// maximum absolute change in elevation + { +#ifdef DEPOS + etrue = efalse = dtrue = dfalse = oppos = 0; +#endif + t += dt; + at = pow(a,theta); + computeFluxes(dt); + double maxc = 0.; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->m = 0; + +#ifdef REDREC +#ifdef DOWN + for ( vector::reverse_iterator p = seq.rbegin(); p != seq.rend(); ++p ) +#else + for ( vector::iterator p = seq.begin(); p != seq.end(); ++p ) +#endif + { + double e = erode ( dt, *p ); + if ( e > maxc) maxc = e; + } +#else + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + double e = erode ( dt, Point(i,j) ); + if ( e > maxc) maxc = e; + } +#endif + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->m = 0; +#ifdef TRANSPORT +#ifdef ICE + removeIceLayer(); +#endif +#ifdef LAKES +#ifndef LAKEFLOWDIR + for ( int i = 0; i < lakes.size(); ++i ) + lakes[i].computeFlowDirection(); + computeWaterLevel(); +#endif +#endif +#ifdef LAYERS + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->adjustLayers(getNode(i,j)->qp-getNode(i,j)->u*dt,t); +#endif +#endif + return maxc; + } + + private: + void adjustWaterLevel ( Point p ) + { + double e, maxc = 0.; + + if ( getNode(p)->m ) return; + getNode(p)->m = 1; + getNode(p)->l = getNode(p)->h; + if ( !getNode(p)->b ) + { + Point dest = getNode(p)->d; + if ( dest == p ) + { + printf ( "Error: Detected node without valid flow target: %i %i\n", + p.i, p.j ); + exit ( -1 ); + } +#ifndef PRECIP + if ( getNode(p)->q == 0 ) + { + printf ( "Error: Detected node without valid flow target: %i %i\n", + p.i, p.j ); + exit ( -1 ); + } +#endif + adjustWaterLevel(dest); + if ( getNode(p)->l < getNode(dest)->l ) getNode(p)->l = getNode(dest)->l; + } + } + + public: + double computeLakeVolume() +// Computes to total volume of all filled lakes. + { + double v = 0.; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) v += getNode(i,j)->l - getNode(i,j)->h; + return v; + } + + void sort ( vector &p ) + { + int i, j, l = p.size()>>1, ir = p.size()-1; + double h; + Point pxt; + + if ( ir == 0 ) return; + for ( ; ; ) + { + if ( l ) h = getNode(pxt=p[--l])->l; + else + { + h = getNode(pxt=p[ir])->l; + p[ir] = p[0]; + if ( !(--ir) ) + { + p[0] = pxt; + return; + } + } + j = ( (i=l) << 1 ) + 1; + while ( j <= ir ) + { + if ( j < ir && getNode(p[j])->l > getNode(p[j+1])->l ) j++; + if ( h > getNode(p[j])->l ) + { + p[i] = p[j]; + j += (i=j)+1; + } + else j = ir+1; + } + p[i] = pxt; + } + } + + template + void sort ( vector> &p ) + { + int i, j, l = p.size()>>1, ir = p.size()-1; + valuetype h; + PointValue pxt; + + if ( ir == 0 ) return; + for ( ; ; ) + { + if ( l ) h = (pxt=p[--l]).d; + else + { + h = (pxt=p[ir]).d; + p[ir] = p[0]; + if ( !(--ir) ) + { + p[0] = pxt; + return; + } + } + j = ( (i=l) << 1 ) + 1; + while ( j <= ir ) + { + if ( j < ir && p[j].d < p[j+1].d ) j++; + if ( h < p[j].d ) + { + p[i] = p[j]; + j += (i=j)+1; + } + else j = ir+1; + } + p[i] = pxt; + } + } + +#ifdef FLEXURE + void restrict ( vector > &src, vector > &dest, + vector > &tmp ) + { + int m = src.size(); + int n = src[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; j += 2 ) + getNode(i,j/2)->w = src[i][(j+n-1)%n]+2*src[i][j]+src[i][(j+1)%n]; + for ( int i = 0; i < m; i += 2 ) + for ( int j = 0; j < n/2; ++j ) + dest[i/2][j] = (getNode((i+m-1)%m,j)->w+2*getNode(i,j)->w+getNode((i+1)%m,j)->w)/16; + } + + void prolongate ( vector > &src, vector > &dest ) + { + int m = src.size(); + int n = src[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) dest[2*i][2*j] = src[i][j]; + m = 2*m; + n = 2*n; + for ( int i = 0; i < m; i += 2 ) + for ( int j = 1; j < n; j += 2 ) + dest[i][j] = 0.5*(dest[i][(j+n-1)%n]+dest[i][(j+1)%n]); + for ( int i = 1; i < m; i += 2 ) + for ( int j = 0; j < n; ++j ) + dest[i][j] = 0.5*(dest[(i+m-1)%m][j]+dest[(i+1)%m][j]); + } + +#ifdef RIGIDITY + void gaussseidel ( vector > &u, vector > &rhs, + double tau, vector > &rig ) + { + double maxchange; + double den = 1 + tau + 20*rig[0][0]; + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( rhs[i][j] < 1e50 ) + u[i][j] = ( rhs[i][j] + + 4*rig[i][j] * ( u[(i+m-1)%m][j] + u[(i+1)%m][j] + + u[i][(j+n-1)%n] + u[i][(j+1)%n] ) + + rig[(i+m-1)%m][j] * ( 4*u[(i+m-1)%m][j] - u[(i+m-2)%m][j] + - u[(i+m-1)%m][(j+n-1)%n] - u[(i+m-1)%m][(j+1)%n] ) + + rig[(i+1)%m][j] * ( 4*u[(i+1)%m][j] - u[(i+2)%m][j] + - u[(i+1)%m][(j+n-1)%n] - u[(i+1)%m][(j+1)%n] ) + + rig[i][(j+n-1)%n] * ( 4*u[i][(j+n-1)%n] - u[i][(j+n-2)%n] + - u[(i+m-1)%m][(j+n-1)%n] - u[(i+1)%m][(j+n-1)%n] ) + + rig[i][(j+1)%n] * ( 4*u[i][(j+1)%n] - u[i][(j+2)%n] + - u[(i+m-1)%m][(j+1)%n] - u[(i+1)%m][(j+1)%n] ) ) + / ( 1 + tau + 16*rig[i][j] + + ( rig[(i+m-1)%m][j] + rig[(i+1)%m][j] + + rig[i][(j+n-1)%n] + rig[i][(j+1)%n] ) ); + } + + void residuum ( vector > &u, vector > &rhs, + vector > &res, double tau, + vector > &rig ) + { + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->w = rig[i][j] * ( u[(i+m-1)%m][j] + u[(i+1)%m][j] + + u[i][(j+n-1)%n] + u[i][(j+1)%n] + - 4*u[i][j] ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + res[i][j] = rhs[i][j] - getNode((i+m-1)%m,j)->w - getNode((i+1)%m,j)->w + - getNode(i,(j+n-1)%n)->w - getNode(i,(j+1)%n)->w + + 4*getNode(i,j)->w - (1+tau)*u[i][j]; + } + +#else + void gaussseidel ( vector > &u, vector > &rhs, + double tau, double rig ) + { + double maxchange; + double den = 1 + tau + 20*rig; + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( rhs[i][j] < 1e50 ) + u[i][j] = ( rhs[i][j] + + 8*rig * ( u[(i+m-1)%m][j] + u[(i+1)%m][j] + + u[i][(j+n-1)%n] + u[i][(j+1)%n] ) + - rig * ( u[(i+m-2)%m][j] + u[(i+2)%m][j] + + u[i][(j+n-2)%n] + u[i][(j+2)%n] + + 2 * ( u[(i+m-1)%m][(j+n-1)%n] + + u[(i+m-1)%m][(j+1)%n] + + u[(i+1)%m][(j+n-1)%n] + + u[(i+1)%m][(j+1)%n] ) ) ) / den; + } + + void residuum ( vector > &u, vector > &rhs, + vector > &res, double tau, double rig ) + { + double diag = 1 + tau + 20*rig; + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->w = rig * ( u[(i+m-1)%m][j] + u[(i+1)%m][j] + + u[i][(j+n-1)%n] + u[i][(j+1)%n] + - 4*u[i][j] ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + res[i][j] = rhs[i][j] - getNode((i+m-1)%m,j)->w - getNode((i+1)%m,j)->w + - getNode(i,(j+n-1)%n)->w - getNode(i,(j+1)%n)->w + + 4*getNode(i,j)->w - (1+tau)*u[i][j]; + } +#endif + + void setzero ( vector > &u ) + { + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) u[i][j] = 0; + } + + void add ( vector > &u, vector > &v ) + { + int m = u.size(); + int n = u[0].size(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) u[i][j] += v[i][j]; + } + + void cycle ( double tau, int l, int ct, int ng, int ngl ) + { + for ( int i = 0; i < ng; ++i ) + gaussseidel ( w[l], frhs[l], tau, rig[l] ); + if ( l+1 < w.size() ) + { + residuum ( w[l], frhs[l], res[l], tau, rig[l] ); + restrict ( res[l], frhs[l+1], res[l] ); + setzero ( w[l+1] ); + for ( int i = 0; i < ct; ++i ) + cycle ( tau, l+1, ct, ng, ngl ); + prolongate ( w[l+1], res[l] ); + add ( w[l], res[l] ); + } + else ng = ngl; + for ( int i = 0; i < ng; ++i ) + gaussseidel ( w[l], frhs[l], tau, rig[l] ); + } + + void applyFlexure ( double dt = 1e99, int ct = 2, int ng = 1, int ngl = 100, int iter = 1 ) + { + double tau = this->tau/dt; + + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *p = getNode(i,j); + p->h -= p->w; + w[0][i][j] = p->w; + frhs[0][i][j] = //p->b ? 1e99 : + -rhoc*(p->h+p->g) + tau*p->w; + } +#ifdef RIGIDITY + if ( rig[0][0][0] != pow(getNode(0,0)->alpha,4)/4 ) + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) rig[0][i][j] = pow(getNode(i,j)->alpha,4)/4; + for ( int l = 1; l < rig.size(); ++l ) + { + restrict ( rig[l-1], rig[l], rig[l-1] ); + for ( int i = 0; i < rig[l].size(); ++i ) + for ( int j = 0; j < rig[l][i].size(); ++j ) rig[l][i][j] /= 16; + } + } +#else + rig[0] = pow(alpha,4)/4; + for ( int l = 1; l < rig.size(); ++l ) rig[l] = rig[l-1]/16; +#endif + for ( int it = 0; it < iter; ++it ) + cycle ( tau, 0, ct, ng, ngl ); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *p = getNode(i,j); + p->w = w[0][i][j]; + p->h += p->w; + } + computeWaterLevel(); + } +#endif + +#ifdef PRECIP + void solveRowCol ( int mode, int n, int i ) +// mode = 0 -> row i; mode = 1 -> col i + { + double beta = (1-lc/ll)*(ll/lf-1); + Matrix dmat(-ld,0,0,-ld); + for ( int j = 0; j < n; ++j ) + { + Node *p = mode ? getNode(j,i) : getNode(i,j); + double f = exp(-p->l/refheight); + Matrix dmat(-ld,0,0,-ld); + if ( !p->b ) + { + if ( refslope > 0 ) + { + double s = p->l - getNode(p->d)->l; + if ( p->d.isDiag(i,j) ) s *= sqrt(0.5); + psi[j] = s > 0. ? exp(-refslope/s) : 0.; + psi[j] = 1-f*evap*(1-psi[j]); + } + else psi[j] = 1-f*evap; + f *= beta/lc; + diag[j] = Matrix(2*ld+1+1/lc,-f-(1-psi[j])/lf,-1/lc,2*ld+1+f+1/lf); + } + else + { + diag[j] = Matrix(1,0,0,1); + dmat = Matrix(0,0,0,0); + f *= lf/(ll-lf); + rhs[j] = Vect(qin/(1+f),qin*f/(1+f)); + } + if ( j > 0 ) + lower[j-1] = dmat; + else + right[0] = dmat; + if ( j < n-1 ) + upper[j] = dmat; + else + bottom[0] = dmat; + } +// for ( int j = 0; j < n-1; ++j ) lower[j] = upper[j] = Matrix(-ld,0,0,-ld); +// right[0] = bottom[0] = Matrix(-ld,0,0,-ld); + for ( int j = 1; j < n-2; ++j ) right[j] = bottom[j] = Matrix(0,0,0,0); + + for ( int j = 0; j < n-1; ++j ) + { + Matrix f = diag[j].inv()*lower[j]; + diag[j+1] -= f*upper[j]; + rhs[j+1] -= f*rhs[j]; + + if ( j < n-2 ) + { + Matrix g = diag[j].inv()*bottom[j]; + if ( j < n-3) + { + right[j+1] -= f*right[j]; + bottom[j+1] -= g*upper[j]; + } + else + { + upper[j+1] -= f*right[j]; + lower[j+1] -= g*upper[j]; + } + diag[n-1] -= g*right[j]; + rhs[n-1] -= g*rhs[j]; + } + } + rhs[n-1] = diag[n-1].inv()*rhs[n-1]; + rhs[n-2] -= upper[n-2]*rhs[n-1]; + rhs[n-2] = diag[n-2].inv()*rhs[n-2]; + for ( int j = n-3; j >= 0; --j ) + { + rhs[j] -= upper[j]*rhs[j+1]; + rhs[j] -= right[j]*rhs[n-1]; + rhs[j] = diag[j].inv()*rhs[j]; + } + for ( int j = 0; j < n; ++j ) + { + if ( rhs[j].u1 < 0 || rhs[j].u2 < 0 ) + { + printf ( "At least one of the fluxes is negative.\n" + "So far I thought that this cannot happen.\n" + "i = %i, j = %i, qv = %e, qc = %e\n", + i, j, rhs[j].u1, rhs[j].u2 ); + exit(1); + } + Node *p = mode ? getNode(j,i) : getNode(i,j); + p->p = psi[j]*rhs[j].u2/lf; +#ifdef EXTPRECIP + p->qv = rhs[j].u1; + p->qc = rhs[j].u2; + p->ptot = rhs[j].u2/lf; +#endif + } + } + + void computePrecipitation() + { + int m = dir < 2 ? this->m : this->n; + int n = dir < 2 ? this->n : this->m; + for ( int i = 0; i < n; ++i ) + { + Node *p = dir < 2 ? getNode ( dir ? m-1 : 0, i ) : getNode ( i, dir > 2 ? m-1 : 0 ); + double f = lf/(ll-lf)*exp(p->l/refheight); + rhs[i] = Vect(qin/(1+f),qin*f/(1+f)); +#ifdef EXTPRECIP + p->qv = rhs[i].u1; + p->qc = rhs[i].u2; +#endif + } + if ( dir%2 ) + for ( int i = m-2; i >= 0; --i ) solveRowCol ( dir > 2, n, i ); + else + for ( int i = 1; i < m; ++i ) solveRowCol ( dir > 1, n, i ); + } +#endif // PRECIP + + void write ( string filename, string mask, int writeheaders = 1, int layer = -1 ) +// Writes the lattice to a file. The mask defines which variables are written. +// For writeheaders = 1, all non-gridded variables are automatically included. +// For writeheaders = 0, only the raw grid data are written. + { + for ( int i = 0; i < mask.size(); ++i ) mask[i] = toupper(mask[i]); + FILE *fp = fopen ( filename.c_str(), "wb" ); + if ( writeheaders ) + { + fwrite ( &m, sizeof(short int), 1, fp ); + fwrite ( &n, sizeof(short int), 1, fp ); + } + map::reverse_iterator it; + for ( it = km.rbegin(); it != km.rend(); ++it ) + { + Keymapentry e = it->second; + if ( e.grid ) + { + int pos = mask.find(e.keystring); + if ( pos != std::string::npos ) + { + for ( int i = pos; i < pos+strlen(e.keystring); ++i ) mask[i] = ' '; +// printf ( "%s\n", mask.c_str() ); + if ( writeheaders ) + { + if ( e.array > 1 && layer < 0 ) + { + fwrite ( "ARR", sizeof(int), 1, fp ); + fwrite ( &e.array, sizeof(int), 1, fp ); + } + fwrite ( &e.keyint, sizeof(int), 1, fp ); + } + if ( strcmp(e.keystring,"S") ) + if ( e.size == 8 && e.fsize == 4 ) + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + float buf = *(double *)((char *)&getNode(i,j)->h+e.offset); + if ( layer < 0 ) + fwrite ( &buf, e.fsize, e.array, fp ); + else + fwrite ( &buf+layer, e.fsize, 1, fp ); + } + else + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( layer < 0 ) + fwrite ( (char *)&getNode(i,j)->h+e.offset, e.fsize, e.array, fp ); + else + fwrite ( (char *)&getNode(i,j)->h+e.offset+layer*e.fsize, e.fsize, 1, fp ); + else + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + float s = getNode(i,j)->l - getNode(getNode(i,j)->d)->l; + if ( getNode(i,j)->d.isDiag(i,j) ) s *= sqrt(0.5); + fwrite ( &s, 4, 1, fp ); + } + } + } + else if ( writeheaders && !e.readonly ) + { + if ( e.array > 1 ) + { + fwrite ( "ARR", sizeof(int), 1, fp ); + fwrite ( &e.array, sizeof(int), 1, fp ); + } + fwrite ( &e.keyint, sizeof(int), 1, fp ); + fwrite ( (char *)&t+e.offset, e.size, e.array, fp ); + } + } + fclose(fp); + } + + void read ( FILE *fp, Keymapentry e, int arr = 1 ) + { + if ( e.size == 8 && e.fsize == 4 ) + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + for ( int k = 0; k < arr; ++k ) + { + float buf; + fread ( &buf, 4, arr, fp ); + double bufd = buf; + if ( k < e.array ) + memcpy ( (char *)&getNode(i,j)->h+e.offset+8*k, &bufd, 8 ); + } + else if ( e.size > 0 ) + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( arr > e.array ) + { + fread ( (char *)&getNode(i,j)->h+e.offset, e.fsize, e.array, fp ); + fseek ( fp, e.fsize*(arr-e.array), SEEK_CUR ); + } + else + fread ( (char *)&getNode(i,j)->h+e.offset, e.fsize, arr, fp ); + else + fseek ( fp, m*n*e.fsize*arr, SEEK_CUR ); + } + + void read ( string filename, string mask, int arr = 0 ) + { + Keymapentry e = getKey(mask.c_str()); + if ( e.grid ) + { + FILE *fp = fopen ( filename.c_str(), "rb" ); + read ( fp, e, arr > 0 ? arr : e.array ); + fclose(fp); + } + else + printf ( "The given key %s could not be assigned to any variable.\n", mask.c_str() ); + } + + double read ( string filename ) +// Reads data from a file and returns time. + { + short int m, n; + int buf; + unsigned int arr = 1, unrec = 0; + FILE *fp = fopen ( filename.c_str(), "rb" ); + fread ( &m, sizeof(short int), 1, fp ); + fread ( &n, sizeof(short int), 1, fp ); + resize(m,n); +#ifdef DX + fread ( dx, sizeof(double), m, fp ); +#endif + while ( fread ( &buf, sizeof(int), 1, fp ) == 1 ) + { + if ( km.find(buf) == km.end() ) + { + char cbuf[5]; + memcpy(cbuf,&buf,4); + cbuf[4] = '\0'; + if ( strcmp(cbuf,"ARR") ) + { + if ( !unrec ) + { + unrec = ftell(fp)-4; + printf ( "Unrecognized key %i at file position %i, looks like '%s'.\n", + buf, unrec, cbuf ); + } + } + else fread ( &arr, sizeof(int), 1, fp ); + } + else + { + Keymapentry e = km[buf]; + if ( unrec ) + { + unrec = 0; + printf ( "Continued reading at file position %i with key '%s'.\n", + ftell(fp)-4, e.keystring ); + } + if ( e.grid ) read(fp,e,arr); + else if ( arr > e.array ) + { + fread ( (char *)&t+e.offset, e.size, e.array, fp ); + fseek ( fp, e.size*(arr-e.array), SEEK_CUR ); + } + else + fread ( (char *)&t+e.offset, e.size, arr, fp ); + arr = 1; + } + } + fclose ( fp ); + return t; + } + + void writeMatlabClass ( string classname = "", string path = "." ) + { + int arrkey; + memcpy(&arrkey,"ARR",4); + if ( classname.length() == 0 ) classname = "openlem"; + FILE *fp = fopen ( (path+"/"+classname+".m").c_str(), "w" ); + fprintf ( fp, "classdef openlem < handle\n\n properties\n m;\n n;\n" ); + char varname[32]; + string varlist; + map::reverse_iterator it; + for ( it = km.rbegin(); it != km.rend(); ++it ) + { + Keymapentry e = it->second; + if ( strcmp(e.keystring,"D") ) + sprintf(varname," %s;\n",e.varname); + else + sprintf(varname," di;\n dj;\n"); + if ( varlist.find(varname) == string::npos ) + varlist.append(varname); + } + fprintf ( fp, "%s", varlist.c_str() ); + fprintf ( fp, " end\n\n methods\n" + " function o = openlem(filename)\n" + " fid = fopen(filename);\n" + " if fid < 0\n" + " o.m = 0;\n" + " o.n = 0;\n" + " return\n" + " end\n" + " arr = 1;\n" + " o.m = fread(fid,1,'int16=>int32');\n" + " o.n = fread(fid,1,'int16=>int32');\n" + " while true\n" + " if arr > 1\n" + " arrsize = [arr,o.n,o.m];\n" + " perm = [ 3 2 1 ];\n" + " else\n" + " arrsize = [o.n,o.m];\n" + " perm = [ 2 1 ];\n" + " end\n" + " key = fread(fid,1,'uint32=>uint32');\n" + " if numel(key) == 0\n" + " break\n" + " end\n" + " switch key\n" ); + for ( it = km.rbegin(); it != km.rend(); ++it ) + { + Keymapentry e = it->second; + fprintf ( fp, " case %u\n", e.keyint ); + if ( !strcmp(e.keystring,"D" ) ) + { + fprintf ( fp, " tmp = fread(fid,[2,o.m*o.n],'uint16')+1;\n" + " o.di = reshape(tmp(1,:),[o.n,o.m])';\n" + " o.dj = reshape(tmp(2,:),[o.n,o.m])';\n" ); + continue; + } +#ifndef PRECIP + if ( !strcmp(e.keystring,"Q" ) ) + { + fprintf ( fp, " o.q = fread(fid,[o.n,o.m],'uint32=>uint32')';\n" ); + continue; + } +#endif + if ( e.grid ) + fprintf ( fp, " o.%s = permute(reshape(fread(fid,prod(arrsize),'", e.varname ); + else + fprintf ( fp, " o.%s = fread(fid,arr,'", e.varname ); + switch ( e.fsize ) + { + case 1: + fprintf ( fp, "int8" ); + break; + case 4: + fprintf ( fp, e.grid ? "float" : "uint32" ); + break; + case 8: + fprintf ( fp, "double" ); + break; + } + fprintf ( fp, "=>" ); + switch ( e.size ) + { + case 1: + fprintf ( fp, "int8" ); + break; + case 8: + fprintf ( fp, "double" ); + break; + default: + fprintf ( fp, e.grid ? "float" : "uint32" ); + } + fprintf ( fp, "'" ); + if ( e.grid ) fprintf ( fp, "),arrsize),perm" ); + fprintf ( fp, ");\n" + " arr = 1;\n" ); + } + fprintf ( fp, " case %u\n", arrkey ); + fprintf ( fp, " arr = fread(fid,1,'uint32');\n" + " end\n" + " end\n" + " fclose(fid);\n" + " end\n\n" + " function write(o,filename)\n" + " fid = fopen(filename,'w');\n" + " fwrite(fid,o.m,'uint16');\n" + " fwrite(fid,o.n,'uint16');\n" ); + for ( it = km.rbegin(); it != km.rend(); ++it ) + { + Keymapentry e = it->second; + if ( e.readonly || strchr(e.keystring,'4') != NULL || strchr(e.keystring,'8') ) + continue; + if ( !strcmp(e.keystring,"D" ) ) + { + fprintf ( fp, " if numel(o.di)\n" ); + fprintf ( fp, " fwrite(fid,%u,'uint32');\n", e.keyint ); + fprintf ( fp, " tmp = [ reshape(o.di',[1,o.m*o.n]); reshape(o.dj',[1,o.m*o.n]) ]-1;\n" + " fwrite(fid,tmp,'uint16')';\n" + " end\n" ); + continue; + } + fprintf ( fp, " if numel(o.%s)\n", e.varname ); + if ( e.grid ) + fprintf ( fp, " arr = 1;\n" + " if numel(size(o.%s)) == 3\n" + " arr = size(o.%s,3);\n" + " end\n", e.varname, e.varname ); + else + fprintf ( fp, " arr = numel(o.%s);\n", e.varname ); + fprintf ( fp, " if arr > 1\n" + " fwrite(fid,%u,'uint32');\n" + " fwrite(fid,arr,'uint32');\n" + " end\n", arrkey ); + fprintf ( fp, " fwrite(fid,%u,'uint32');\n", e.keyint ); +#ifndef PRECIP + if ( !strcmp(e.keystring,"Q" ) ) + { + fprintf ( fp, " fwrite(fid,o.q','uint32')';\n" + " end\n" ); + continue; + } +#endif +// fprintf ( fp, " fwrite(fid,o.%s','", e.varname ); + if ( e.grid ) + fprintf ( fp, " fwrite(fid,permute(o.%s,numel(size(o.%s)):-1:1),'" , e.varname, e.varname ); + else + fprintf ( fp, " fwrite(fid,o.%s,'", e.varname ); + switch ( e.fsize ) + { + case 1: + fprintf ( fp, "int8" ); + break; + case 4: + fprintf ( fp, e.grid ? "float" : "uint32" ); + break; + case 8: + fprintf ( fp, "double" ); + break; + } + fprintf ( fp, "');\n" + " end\n" ); + } + fprintf ( fp, " fclose(fid);\n" + " end\n\n" + " end\n\n" + "end\n" ); + fclose ( fp ); + } + + class Delta + { + public: + Point first; + Point last; + int i; + double q; + double d; + Grid *g; + + Delta ( Grid *g, Point p ) + { + this->g = g; + first = last = p; + i = 0; + q = g->getNode(p)->q; + d = g->offset[i].d/q; + } + + void add ( Point p ) + { + g->getNode(last)->d = p; + } + + Point point() + { + return Point((first.i+g->offset[i].p.i)%g->m,(first.j+g->offset[i].p.j)%g->n); + } + + void step() + { + d = g->offset[++i].d/q; + } + }; + + class Lake + { + public: + double l; + vector points; + vector boundary; + vector delta; + map> rim; + Point outlet; + Grid *g; + + Node *getNode ( Point p ) + { + return g->getNode(p); + } + + Node *getNode ( vector::iterator p ) + { + return g->getNode(p); + } + + int add ( Point p ) + { + Node *pn = getNode(p); + if ( pn->b || pn->m ) return 0; + pn->m = 1; + points.push_back(p); + int n = 1; + vector neigh = g->getNeighbors(p); + for ( vector::iterator it = neigh.begin(); it != neigh.end(); ++it ) + { + pn = getNode(it); + if ( !pn->m ) + rim[pn->l].insert(*(unsigned int *)&*it); + } + return n; + } + + void addDonors ( Point p ) + { +// printf ( "Add %i %i\n", p.i, p.j ); + vector neigh = g->getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( g->getNode(u)->drainsTo(p) ) + { +// printf ( "Neigh %i %i %f %f\n", u->i, u->j, g->getNode(u)->l, l ); + if ( g->getNode(u)->l == l ) + { + points.push_back(*u); + addDonors(*u); + } + else + boundary.push_back(*u); + } + } + + public: + Lake ( Grid *g ) + { + this->g = g; + } + + Lake ( Grid *g, Point p ) + { + this->g = g; + int lake = 0; + double q = 0; + Node *pn = getNode(p); + l = pn->l; + while ( pn->l == l ) + { + outlet = p; + if ( pn->b ) break; + pn = getNode(p=pn->d); + } + vector neigh = g->getNeighbors(outlet); + for ( vector::iterator p = neigh.begin(); p != neigh.end(); ++p ) + if ( (pn=getNode(p))->drainsTo(outlet) ) + if ( pn->l == l ) + { + points.push_back(*p); + if ( pn->q > q ) + { + q = pn->q; + lake = 1; + } + addDonors(*p); + } + else if ( pn->q > q ) + { + q = pn->q; + lake = 0; + } + double depth = 0; + for ( vector::iterator p = points.begin(); p != points.end(); ++p ) + if ( getNode(p)->l-getNode(p)->h >= depth ) depth = getNode(p)->l-getNode(p)->h; + if ( lake ) + { + for ( vector::iterator p = points.begin(); p != points.end(); ++p ) + { + getNode(p)->l = getNode(p)->q = 0; + getNode(p)->m = 2; + } + int nb = 0; + for ( vector::iterator p = boundary.begin(); p != boundary.end(); ++p ) + { + double h = (pn=getNode(p))->h, maxs = 0.; + vector neigh = g->getNeighbors(*p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(u)->m == 2 ) + { + double s = h-getNode(u)->h; + if ( p->isDiag(*u) ) s *= sqrt(0.5); + if ( s > maxs ) + { + pn->d = *u; + maxs = s; + } + } + if ( getNode(pn->d)->q == 0 ) boundary[nb++] = pn->d; + getNode(pn->d)->q += pn->q; + } + boundary.resize(nb); +// Now boundary should consist of all points that receive discharge from outside. +// new part +// for ( int i = 0; i < nb; ++i ) delta.push_back(Delta(g,boundary[i])); + +// Find biggest delta +// int ilowest = 0; +// for ( int i = 0; i < delta.size(); ++i ) +// if ( delta[i].d < delta[ilowest].d ) ilowest = i; +// + +// end of new part + for ( vector::iterator p = points.begin(); p != points.end(); ++p ) + { + pn = getNode(p); + for ( vector::iterator b = boundary.begin(); b != boundary.end(); ++b ) + pn->l += getNode(b)->q/(g->computeDistSquare(*p,*b)+1.); + } + g->sort(points); + points.push_back(outlet); + for ( vector::iterator p = points.begin()+1; p != points.end(); ++p ) + { + getNode(p-1)->d = *p; + if ( g->computeDistSquare(*(p-1),*p) > 2 ) + getNode(p)->trunk = *(p-1); +// printf ( "(%i,%i) -> (%i,%i) %i (%i,%i)\n", (p-1)->i, (p-1)->j, p->i, p->j, g->computeDistSquare(*(p-1),*p), getNode(p)->trunk.i, getNode(p)->trunk.j ); + } +// printf ( "%i %i %f %f\n", p->j+1, p->i+1, g->getNode(p)->q, g->getNode(p)->l ); +// getNode(points.back())->d = outlet; + for ( vector::iterator p = points.begin(); p != points.end(); ++p ) + { + getNode(p)->l = l; + getNode(p)->m = 3; + } + printf ( "Lake " ); + } + else + { + for ( vector::iterator p = points.begin(); p != points.end(); ++p ) + { + getNode(p)->l = getNode(p)->h; + getNode(p)->m = 1; + } + printf ( "Floodplain " ); + } + printf ( "outlet = (%i,%i), area = %i, depth = %f\n", outlet.i, outlet.j, points.size(), depth ); + } + + /* + void printRim() + { + for ( map>::iterator it = rim.begin(); it != rim.end(); ++it ) + for ( set::iterator q = it->second.begin(); q != it->second.end(); ++q ) + { + Point xp = *(Point*)&*q; + printf ( "%e -> (%i,%i) %i\n", it->first, xp.i, xp.j, *q ); + } + } + */ + + int fill ( Point p ) + { + int n; + + points.clear(); + rim.clear(); + if ( add(p) == 0 ) return 0; + n = 1; + l = getNode(p)->l; + set lowest; + while ( 1 ) + { + double minelev = rim.begin()->first; + lowest = rim.begin()->second; + if ( minelev < l ) break; + rim.erase(rim.begin()); + for ( set::iterator q = lowest.begin(); q != lowest.end(); ++q ) + n += add(*(Point *)&*q); + if ( minelev == l ) + while ( rim.count(l) ) + { + lowest = rim[l]; + rim.erase(l); + for ( set::iterator q = lowest.begin(); q != lowest.end(); ++q ) + n += add(*(Point *)&*q); + } + l = minelev; + } + outlet = *(Point *)&*(lowest.begin()); + for ( int i = 0; i < points.size(); ++i ) + getNode(points[i])->l = l; + return points.size(); + } + + void computeFlowDirection() + { + Node *p; + int osize = points.size(); + for ( int i = 0; i < points.size(); ++i ) getNode(points[i])->m = 1; + points.clear(); + points.push_back(outlet); + getNode(outlet)->m = 0; + boundary.clear(); + for ( int i = 0; i < points.size(); ++i ) + { + vector neigh = g->getNeighbors(points[i]); + for ( vector::iterator p = neigh.begin(); p != neigh.end(); ++p ) + if ( getNode(p)->m ) + { + getNode(p)->m = 0; + if ( !getNode(p)->b) getNode(p)->d = points[i]; + getNode(p)->l = getNode(p)->h; + if ( getNode(p)->l <= getNode(points[i])->l ) + { + getNode(p)->l = getNode(points[i])->l; + points.push_back(*p); + } + else + { +// printf ( "Silted up: (%i,%i)\n", p->i, p->j ); + boundary.push_back(*p); + } + } + if ( i == points.size()-1 && boundary.size() > 0 ) + { + int jmin = 0; + for ( int j = 1; j < boundary.size(); ++j ) + if ( getNode(boundary[j])->h < getNode(boundary[jmin])->h ) jmin = j; + points.push_back(boundary[jmin]); +// printf ( "Added silted-up point (%i,%i) %i\n", boundary[jmin].i, boundary[jmin].j, points.size() ); + boundary[jmin] = boundary.back(); + boundary.pop_back(); + } + } + for ( int i = 0; i < points.size(); ++i ) getNode(points[i])->m = 0; + if ( points.size() < osize ) + printf ( "Lost lake points %i -> %i\n", osize, points.size() ); + } + + void clear() + { + points.clear(); + boundary.clear(); + } + }; + + void fillLakes() +// Fills all local depressions by determining a local water level. +// Flow directions of all points are computed. + { + Lake lake(this); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->l = getNode(i,j)->h; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + getNode(i,j)->d.set(i,j); + if ( getNode(i,j)->b == 0 ) + computeFlowDirection( Point(i,j) ); + } + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->b == 0 && getNode(i,j)->d.equals(i,j) ) + { + lake.fill(Point(i,j)); + lake.computeFlowDirection(); + lake.clear(); + } + printf ( "Finished\n" ); + computeFlowDirection(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->m ) + { + printf ( "Unresolved marker: %i %i\n", i, j ); + exit ( 0 ); + } + } + + void findLakes() + { + Node *pn; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + getNode(i,j)->trunk = Point(i,j); +#ifndef ICE + computeDischarge(); +// Otherwise it has already been computed +#endif + printf ( "WL2 %i\n", computeWaterLevel() ); + lakes.clear(); + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( !(pn=getNode(i,j))->m && !pn->b && pn->l-pn->h > mld ) + { + Lake tmp(this,Point(i,j)); + if ( getNode(tmp.outlet)->m == 3 ) lakes.push_back(tmp); + } + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + if ( (pn=getNode(i,j))->m != 3 ) pn->l = pn->h; + pn->m = 0; + } +// for ( int i = 0; i < lakes.size(); ++i ) +// printf ( "Lake %i, outlet = (%i,%i), area = %i\n", i, lakes[i].outlet.i, lakes[i].outlet.j, lakes[i].points.size() ); + } + + }; + + class Delta + { + public: + Point first; + Point last; + int i; + double q; + double d; + Grid *g; + + Delta ( Grid *g, Point p ) + { + this->g = g; + first = last = p; + i = 0; + q = g->getNode(p)->q; + d = g->offset[i].d/q; + } + + void add ( Point p ) + { + g->getNode(last)->d = p; + } + + Point point() + { + return Point((first.i+g->offset[i].p.i)%g->m,(first.j+g->offset[i].p.j)%g->n); + } + + void step() + { + d = g->offset[++i].d/q; + } + }; +} diff --git a/include/aspect/mesh_deformation/openlem.h b/include/aspect/mesh_deformation/openlem.h new file mode 100644 index 00000000000..c1176597357 --- /dev/null +++ b/include/aspect/mesh_deformation/openlem.h @@ -0,0 +1,900 @@ +/* + Copyright (C) 2018 - 2024 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . +*/ + + +#ifndef _aspect_mesh_deformation_openlem_h +#define _aspect_mesh_deformation_openlem_h + +#include +#include +#include +#include +#include +#include + +//#include +#include + +namespace openlem +{ + class OceanGrid : public Grid + { + public: + double odiff; + + OceanGrid ( int m = 1, int n = 1 ) : Grid(m,n) + { + addKey ( "od", "od", &odiff, sizeof(odiff), 8, 0 ); + } + + void clearOceans() + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) getNode(i,j)->b &= 1; + } + + void markOcean ( Point p, double l ) + { + if ( getNode(p)->h <= l ) + { + getNode(p)->l = l; + getNode(p)->b |= 2; + vector neigh = getNeighbors(p); + for ( vector::iterator u = neigh.begin(); u != neigh.end(); ++u ) + if ( getNode(u)->drainsTo(p) ) + markOcean(*u,l); + } + } + + vector > findDeltas ( double dt ) + { + return emplaceSediments(dt); + } + + vector > emplaceSediments ( double dt ) + { + vector > delta; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->b & 2 ) + { + if ( getNode(i,j)->qs ) + delta.push_back(PointValue(Point(i,j),getNode(i,j)->qs)); + getNode(i,j)->h += getNode(i,j)->u*dt; + } + sort(delta); + if ( offset.size() < m*n ) + { + printf ( "Computing offsets\n" ); + offset.resize(m*n); + int k = 0; + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + offset[k].p = Point(i,j); + int di = i <= m/2 ? i : m-i; + int dj = j <= n/2 ? j : n-j; + offset[k++].d = di*di+dj*dj; + } + offset[0].d = -1; + sort(offset); + } + for ( int i = 0; i < delta.size(); ++i ) + { + double v = delta[i].d*dt; + double tmp; + double vd = 0; + Point p = delta[i].p; + vector >::iterator it = offset.begin(); + while ( it != offset.end() ) + { + if ( (tmp=getNode(p)->h+v) <= getNode(p)->l ) + { +#ifdef LAYERS + getNode(p)->adjustLayers(v); +#endif + getNode(p)->h = tmp; + vd += v; + break; + } + else + { + v = v-(tmp=getNode(p)->l-getNode(p)->h); +#ifdef LAYERS + getNode(p)->adjustLayers(tmp); +#endif + getNode(p)->h = getNode(p)->l; + do + { + p = (*(it++)).p; + p = Point((delta[i].p.i+p.i)%m,(delta[i].p.j+p.j)%n); + } + while ( (getNode(p)->b&2) == 0 ); + } + } + } + return delta; + } + +#ifdef LAYERS + void diffuse ( double dt, Point p ) + { + Node *pn = getNode(p); + double dhsum = 0; + vector neigh = getNearestNeighbors(p); + vector dh(4); + for ( int i = 0; i < 4; ++i ) + if ( getNode(neigh[i])->m ) + if ( ( dh[i] = pn->qp-getNode(neigh[i])->qp ) < 0 ) + { + dh[i] = 0; + diffuse ( dt, neigh[i] ); + } + else + dhsum += dh[i]; + dhsum *= dt*odiff; + double avail = pn->h-pn->qp; + if ( pn->bottom[0] < 0) avail -= pn->bottom[0]; + double f = 1; + if ( dhsum > avail ) f = avail/dhsum; + for ( int i = 0; i < 4; ++i ) + getNode(neigh[i])->h += f*dt*odiff*dh[i]; + pn->h -= f*dhsum; + pn->adjustLayers(pn->h-pn->qp,t); + pn->m = 0; + } + + void diffuse ( double dt ) + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + if ( pn->b&2 ) pn->qp = pn->h; + pn->m = pn->b&2; + } + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + if ( getNode(i,j)->m ) diffuse(dt,Point(i,j)); + } + +#else + + void diffuse ( double dt ) + { + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + if ( pn->b&2 ) pn->qp = pn->h; + } + for ( int i = 0; i < m; ++i ) + for ( int j = 0; j < n; ++j ) + { + Node *pn = getNode(i,j); + if ( pn->b&2 ) + { + if ( getNodeP(i+1,j)->b&2 ) + { + double dh = dt*odiff*(pn->qp-getNodeP(i+1,j)->qp); + pn->h -= dh; + getNodeP(i+1,j)->h += dh; + } + if ( getNodeP(i,j+1)->b&2 ) + { + double dh = dt*odiff*(pn->qp-getNodeP(i,j+1)->qp); + pn->h -= dh; + getNodeP(i,j+1)->h += dh; + } + } + } + } +#endif + + + }; + + class Connector + { + public: + double hscale, x0, y0, alpha, dx0, dy0, dalpha; + vector> x, y, vx, vy, vz; + OceanGrid *g; + + Connector ( OceanGrid *g, double hscale = 1, double x0 = 0, double y0 = 0, double alpha = 0 ) + { + std::cout << "connector construct << hscale = " << hscale << ", x0:y0 = " << x0 << ":" << y0 << ", alpha = " << alpha << std::endl; + this->g = g; + this->hscale = hscale; + this->x0 = x0; + this->y0 = y0; + this->alpha = alpha; + dx0 = dy0 = dalpha = 0; + x.resize(g->m); + for ( int i = 0; i < g->m; ++i ) x[i].resize(g->n); + y.resize(g->m); + for ( int i = 0; i < g->m; ++i ) y[i].resize(g->n); + vx.resize(g->m); + for ( int i = 0; i < g->m; ++i ) vx[i].resize(g->n); + vy.resize(g->m); + for ( int i = 0; i < g->m; ++i ) vy[i].resize(g->n); + vz.resize(g->m); + for ( int i = 0; i < g->m; ++i ) vz[i].resize(g->n); + + updateXY(); + } + + void updateXY() + { + std::cout << "connector updateXY" << std::endl; + double c = cos(alpha); + double s = sin(alpha); + for ( int i = 0; i < g->m; ++i ) + for ( int j = 0; j < g->n; ++j ) + { + x[i][j] = x0 + hscale*(c*i-s*j); + y[i][j] = y0 + hscale*(s*i+c*j); + //if (i == 0 && j == 0) + // { + // std::cout << "x:y = " << x[i][j] << ":" << y[i][j] << ", x0:y0 = " << x0 << ":" << y0 << ", hscale = " << hscale << std::endl; + // } + } + } + + void convertVelocities(bool minimize_advection = true) + { + + //todo: add and substract 4 of b + std::cout << "connector convertVelocities" << std::endl; + int n = 0; + double xmean = 0, ymean = 0, vxmean = 0, vymean = 0, cross = 0, rsq = 0; + double c = cos(alpha); + double s = sin(alpha); + + for ( int i = 0; i < g->m; ++i ) + for ( int j = 0; j < g->n; ++j ) + if (!std::isnan(vx[i][j])) + { +// Convert velocities to OpenLEM coordinate system + double tmp = vx[i][j]; + vx[i][j] = (c*vx[i][j]+s*vy[i][j])/hscale;; + vy[i][j] = (-s*tmp+c*vy[i][j])/hscale; + + AssertThrow(!std::isnan(vx[i][j]) && !std::isnan(vy[i][j]), + aspect::ExcMessage("vx and/or vy are nan at " + std::to_string(i) + ":" + std::to_string(j) + ", vxx:y = " + std::to_string(vx[i][j]) + ":" + std::to_string(vy[i][j]))); +// Only continental area + if ( g->getNode(i,j)->b == 0 ) + { + ++n; + xmean += i; + ymean += j; + vxmean += vx[i][j]; + vymean += vy[i][j]; + cross += i*vy[i][j]-j*vx[i][j]; + //if (std::fabs(vx[i][j]-(-y[i][j])) > std::numeric_limits::epsilon() || std::fabs(vy[i][j]-x[i][j]) > std::numeric_limits::epsilon()) + // { + // std::cout << "wrong: i:j = " + std::to_string(i) + ":" + std::to_string(j) + ", x:y = " << x[i][j] << ":" << y[i][j] << ", vx:vy = " < 0); + xmean /= n; + ymean /= n; + vxmean /= n; + vymean /= n; + cross /= n; + rsq /= n; + printf ( "%e %e\n", xmean, ymean ); + dalpha = (cross-xmean*vymean+ymean*vxmean)/(rsq-xmean*xmean-ymean*ymean); + std::cout << "cross = " << cross << ", rsq = " << rsq << ", cross-rsq = " << cross-rsq << ", mean diff top = " << -xmean *vymean+ymean *vxmean << ", mean diff bottom = " << -xmean *xmean-ymean *ymean << std::endl; + dx0 = vxmean+dalpha*ymean; + dy0 = vymean-dalpha*xmean; + printf ( "%e %e %e\n", dalpha, dx0, dy0 ); + //std::cout << " vx:vy = " << 0 << ":" << 0 << " = "<< vx[0][0] << " : " << vy[0][0] << std::endl; + //for ( int i = 0; i < g->m; i=i+10 ) + // for ( int j = 0; j < g->n; j=j+10 ) + // { + // if ( + // //(i == 5 && j == 5) || + // //(i == 50 && j == 50) || + // //(i == 100 && j == 100) || + // (i == 116 && j == 26) + // ) + // { + // std::cout << " vx:vy = " << i << ":" << j << " = "<< vx[i][j] << " : " << vy[i][j] << std::endl; + // } + // } + for ( int i = 0; i < g->m; ++i ) + for ( int j = 0; j < g->n; ++j ) + if (!std::isnan(vx[i][j])) + { + //if ( + // //(i == 5 && j == 5) || + // //(i == 50 && j == 50) || + // //(i == 100 && j == 100) || + // (i == 116 && j == 26) + //) + // { + // std::cout << " i:j= " << i << ":" << j << " vx:vy = "<< vx[i][j] << " : " << vy[i][j] << ", vx-dx0-dalpha*j = " << vx[i][j]-(dx0-dalpha *j) << ", vy-dy0+daplha*i = " << vy[i][j]-(dy0+dalpha *i) << std::endl; + // } + vx[i][j] -= dx0-dalpha*j; + vy[i][j] -= dy0+dalpha*i; + AssertThrow(!std::isnan(vx[i][j]) && !std::isnan(vy[i][j]), + aspect::ExcMessage("vx and/or vy are nan at " + std::to_string(i) + ":" + std::to_string(j) + ", vxx:y = " + std::to_string(vx[i][j]) + ":" + std::to_string(vy[i][j]))); + } + + double dx0_old = dx0; + dx0 = (c*dx0-s*dy0)*hscale; + dy0 = (s*dx0_old+c*dy0)*hscale; + + //std::cout << "residual x:y = " << 0 << ":" << 0 << " = "<< vx[0][0] i << " : " << vy[0][0] << std::endl; + //std::cout << "residual x:y = " << 0 << ":" << openlem_ny << " = "<< vx[0][openlem_ny-1] << " : " << vy[0][openlem_ny-1] << std::endl; + //std::cout << "residual x:y = " << 0 << ":" << 0 << " = "<< vx[openlem_nx-1][openlem_ny-1] << " : " << vy[openlem_nx-1][openlem_ny-1] << std::endl; + //std::cout << "residual x:y = " << 0 << ":" << 0 << " = "<< vx[openlem_nx-1][0] << " : " << vy[openlem_nx-1][0] << std::endl; + //for ( int i = 0; i < g->m; i=i+10 ) + // for ( int j = 0; j < g->n; j=j+10 ) + // { + // //if (vx[i][j] > 0 && vy[i][j] > 0) + // if ( + // //(i == 5 && j == 5) || + // //(i == 50 && j == 50) || + // //(i == 100 && j == 100) || + // (i == 116 && j == 126) + // ) + // { + // std::cout << "residual x:y = " << i << ":" << j << " = "<< vx[i][j] << " : " << vy[i][j] << std::endl; + // } + // } + } + + void computeUpliftRate() + { + std::cout << "connector computeUpliftRate" << std::endl; + for ( int i = 0; i < g->m; ++i ) + for ( int j = 0; j < g->n; ++j ) + if ( g->getNode(i,j)->b &1 && !std::isnan(vx[i][j])) + { + double dhdx = vx[i][j]<0 ? + g->getNodeP(i+1,j)->h-g->getNode(i,j)->h : + g->getNode(i,j)->h-g->getNodeP(i-1,j)->h; + double dhdy = vy[i][j]<0 ? + g->getNodeP(i,j+1)->h-g->getNode(i,j)->h : + g->getNode(i,j)->h-g->getNodeP(i,j-1)->h; + //if (i == 19 && j == 24) + // std::cout << "before dhdx:y = " << dhdx << ':' << dhdy << ", u = " << g->getNode(i,j)->u << ", vx:y:z = " << vx[i][j] << ":" << vy[i][j] << ":" << vz[i][j] << std::endl; + g->getNode(i,j)->u = vz[i][j]-vx[i][j]*dhdx-vy[i][j]*dhdy; + //if (i == 19 && j == 24) + // std::cout << "dhdx:y = " << dhdx << ':' << dhdy << ", u = " << g->getNode(i,j)->u << ", vx:y:z = " << vx[i][j] << ":" << vy[i][j] << ":" << vz[i][j] << std::endl; + } + } + + void updateCoordinateSystem ( double dt ) + { + std::cout << ".. connector updateCoordinatesystem .."; + //std::cout << "connector updateCoordinatesystem x0:y0 = " << x0 << ":" << y0 << ", dx0:dy0 = " << dx0 << ":" << dy0 << ", dt = " << dt << ", aplha = " << alpha << ", dalhpa = " << dalpha <= interperolat.size()) + { + i = interperolat.size()-1; + } + int j = floor(y); + if (j < 0) + { + j = 0; + } + else if (j >= interperolat[i].size()) + { + j = interperolat[i].size()-1; + }*/ + // find the closest point + double closest_distance = std::numeric_limits::infinity(); + size_t closest_xi = std::numeric_limits::signaling_NaN(); + size_t closest_yi = std::numeric_limits::signaling_NaN(); + for (unsigned int xi = 0; xi < interperolat.size()-1; ++xi) + { + for (unsigned int yi = 0; yi < interperolat[xi].size()-1; ++yi) + { + double distance = (x-xi*dx-origin_x)*(x-xi*dx-origin_x)+(y-yi*dy-origin_y)*(y-yi*dy-origin_y); + if (distance < closest_distance) + { + closest_distance = distance; + closest_xi = xi; + closest_yi = yi; + } + + } + } + + //if (x > -6345.87-100. && x < -6345.87+4000. && y > 28103. - 200. && y < 28103. + 200. ) + //if (x > -6345.87-100. && x < -6345.87+5000. && y > 28103. - 400. && y < 28103. + 400. ) + //if (x > 3125-100. && x < 3125+5000. && y > 45312.5 - 400. && y < 45312.5 + 400. ) + //std::cout << "closest_xi = " << closest_xi << ", closest_yi = " << closest_yi << ", interp = " << interperolat[closest_xi][closest_yi] << ", interp2 = " << interperolat2[closest_xi][closest_yi] << std::endl; + return interperolat[closest_xi][closest_yi];// - interperolat2[closest_xi][closest_yi];//g->getNode(i,j)->h;//interperolat[i][j]; + } + + void write_vtk(double reference_surface_height, int timestep, double time, std::string path, std::string prestring = "") const + { + std::string timestep_string = std::to_string(timestep); + size_t n_zero = 5; + std::string relative_path = "openLEM/openlem_surface" + prestring + std::string(n_zero - std::min(n_zero, timestep_string.length()), '0') +timestep_string +".vtu"; + std::string full_path = path + relative_path; + std::string pvtu_full_path = path + "openlem_surface" + prestring + ".pvd"; + std::cout << "writing openlem VTK file to " << full_path << std::endl; + std::vector grid_x(0); + std::vector grid_y(0); + std::vector grid_z(0); + std::vector grid_elevation(0); + std::vector grid_uplift(0); + std::vector grid_boundary(0); + std::vector grid_q(0); + std::vector grid_i(0); + std::vector grid_j(0); + std::vector grid_vx(0); + std::vector grid_vy(0); + + std::vector> grid_connectivity(0); + + unsigned int n_cell = x.size() * y.size(); + unsigned int n_p = (x.size() + 1) * (y.size() + 1); + + + grid_x.resize(n_p); + grid_z.resize(n_p); + grid_y.resize(n_p); + + grid_elevation.resize(n_p); + grid_uplift.resize(n_p); + grid_boundary.resize(n_p); + grid_q.resize(n_p); + grid_i.resize(n_p); + grid_j.resize(n_p); + grid_vx.resize(n_p); + grid_vy.resize(n_p); + grid_connectivity.resize(n_cell,std::vector((2-1)*4)); + + size_t counter = 0; + for (size_t j = 0; j <= y.size()-1; ++j) + { + for (size_t i = 0; i <= x.size()-1; ++i) + { + grid_x[counter] = x[i][j]; + grid_y[counter] = y[i][j]; + grid_z[counter] = reference_surface_height + g->getNode(i,j)->h; + grid_elevation[counter] = g->getNode(i,j)->h; + grid_uplift[counter] = g->getNode(i,j)->u; + grid_boundary[counter] = g->getNode(i,j)->b; + grid_q[counter] = g->getNode(i,j)->q; + grid_i[counter] = i; + grid_j[counter] = j; + grid_vx[counter] = vx[i][j];//g->getNode(i,j)->u; + grid_vy[counter] = vy[i][j];//g->getNode(i,j)->u; + counter++; + } + } + + counter = 0; + for (size_t j = 1; j <= y.size()-1; ++j) + { + for (size_t i = 1; i <= x.size()-1; ++i) + { + grid_connectivity[counter][0] = i + (j - 1) * (x.size()-1 + 1) - 1; + grid_connectivity[counter][1] = i + 1 + (j - 1) * (x.size()-1 + 1) - 1; + grid_connectivity[counter][2] = i + 1 + j * (x.size()-1 + 1) - 1; + grid_connectivity[counter][3] = i + j * (x.size()-1 + 1) - 1; + counter++; + } + } + std::stringstream buffer; + std::ofstream myfile; + myfile.open (full_path); + buffer << " " << std::endl; + buffer << R"()" << std::endl; + buffer << "" << std::endl; + buffer << "" << std::endl; + buffer << R"(0)" << std::endl; + buffer << "" << std::endl; + buffer << "" << std::endl; + buffer << " " << std::endl; + buffer << R"( )" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_x[i] << " " << grid_y[i] << " " << grid_z[i] << std::endl; + } + buffer << " " << std::endl; + buffer << " " << std::endl; + buffer << std::endl; + buffer << " " << std::endl; + buffer << R"( )" << std::endl; + for (size_t i = 0; i < n_cell; ++i) + buffer << grid_connectivity[i][0] << " " <" << std::endl; + buffer << R"( )" << std::endl; + for (size_t i = 1; i <= n_cell; ++i) + buffer << i * 4 << " "; + buffer << std::endl << " " << std::endl; + buffer << R"( )" << std::endl; + for (size_t i = 0; i < n_cell; ++i) + buffer << "9" << " "; + buffer << std::endl <<" " << std::endl; + buffer << " " << std::endl; + + buffer << " " << std::endl; + + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_elevation[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_uplift[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_boundary[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_q[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_i[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_j[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_vx[i] << std::endl; + } + buffer << "" << std::endl; + buffer << R"()" << std::endl; + for (size_t i = 0; i < n_p; ++i) + { + buffer << grid_vy[i] << std::endl; + } + buffer << "" << std::endl; + buffer << " " << std::endl; + //buffer << " " << std::endl; + //buffer << R"()" << std::endl; + + //for (size_t i = 0; i < n_p; ++i) + // { + // buffer << 1 << " " << 2 << " " << 3 << std::endl; + // } + //buffer << "" << std::endl; + buffer << " " << std::endl; + buffer << " " << std::endl; + buffer << "" << std::endl; + myfile << buffer.str(); + buffer.str(std::string()); + + myfile.close(); + if (prestring != "") + { + std::cout << "Finished writing openlem VTK file" << std::endl; + return; + } + std::string new_pvtu_line = " "; + if (timestep == 0) + { + + myfile.open (pvtu_full_path); + + buffer << R"()" << std::endl; + buffer << R"()"; + buffer << R"( )" << std::endl; + buffer << new_pvtu_line << std::endl; + buffer << R"( )" << std::endl; + buffer << R"()" << std::endl; + myfile << buffer.str(); + buffer.str(std::string()); + myfile.close(); + } + else + { + // read file into vector + std::ifstream infile(pvtu_full_path); + std::string line; + std::vector lines_vector; + while (std::getline(infile, line)) + { + std::istringstream iss(line); + lines_vector.emplace_back(line); + } + infile.close(); + lines_vector.insert(lines_vector.end() - 2, new_pvtu_line + '\n'); + + myfile.open (pvtu_full_path); + + for (std::string &line : lines_vector) + { + buffer << line << std::endl; + } + myfile << buffer.str(); + buffer.str(std::string()); + } + + //std::cout << " \r"; + //std::cout.flush(); + std::cout << "Finished writing openlem VTK file" << std::endl; + } + }; +} +namespace aspect +{ + namespace MeshDeformation + { + /** + * A class that represents a mesh deformation function that can be + * prescribed on the boundary of the domain. + * + * @ingroup MeshDeformation + */ + template + class OpenLEM : public Interface, public SimulatorAccess + { + public: + struct PositionVelocity + { + int x; + int y; + Tensor<1,dim> velocity; + }; + /** + * Constructor. + */ + OpenLEM(); + + /** + * Initialize variables for openlem. + */ + virtual void initialize () override; + + /** + * + */ + void update() override; + + /** + * A function that creates constraints for the velocity of certain mesh + * vertices (e.g. the surface vertices) for a specific boundary. + * The calling class will respect + * these constraints when computing the new vertex positions. + */ + void + compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_id) const override; + + /** + * Returns whether or not the plugin requires surface stabilization + */ + bool needs_surface_stabilization () const override; + + /** + * Declare parameters for the free surface handling. + */ + static + void declare_parameters (ParameterHandler &prm); + + /** + * Parse parameters for the free surface handling. + */ + void parse_parameters (ParameterHandler &prm) override; + + private: + /** + * Execute openlem + */ + void execute_openlem(openlem::OceanGrid &grid, + //std::vector &elevation, + //std::vector &extra_vtk_field, + //std::vector &velocity_x, + //std::vector &velocity_y, + //std::vector &velocity_z, + const double &openlem_timestep_in_years, + const unsigned int &openlem_iterations, + const std::string &dirname); + + /** + * Function to fill the openlem arrays (height and velocities) with the data received from ASPECT in the correct index order. + */ + void fill_openlem_arrays(openlem::OceanGrid &grid, + // std::vector &elevation, + // std::vector &bedrock_transport_coefficient_array, + // std::vector &bedrock_river_incision_rate_array, + // std::vector &velocity_x, + // std::vector &velocity_y, + // std::vector &velocity_z, + std::vector> &temporary_variables); + + /** + * Function to get the ASPECT topography and velocities at the surface, and an index for transferring these to openlem. + */ + std::vector> get_aspect_values() const; + + /** + * define a old and new grid, so that we can compute a difference + */ + openlem::OceanGrid grid_old; + openlem::OceanGrid grid_new; + openlem::Connector connector; + unsigned int openlem_nx; + unsigned int openlem_ny; + unsigned int openlem_dx; + unsigned int openlem_dy; + unsigned int openlem_x_extent; + unsigned int openlem_y_extent; + unsigned int aspect_nx; + unsigned int aspect_ny; + double aspect_dx; + double aspect_dy; + unsigned int aspect_x_extent; + unsigned int aspect_y_extent; + openlem::Point deepest_point; + double openlem_ocean_diffusivity; + + bool openlem_minimize_advection; + double openlem_kd; + double openlem_kt; + //std::vector> mesh_velocity_z; + std::vector> aspect_mesh_dh; + std::vector> aspect_mesh_z; + std::vector> openlem_mesh_aspect_cell; + std::vector> aspect_mesh_aspect_cell_id; + std::vector> aspect_mesh_aspect_cell_face_number; + /** + * Variable to hold ASPECT domain extents. + */ + std::array,dim> grid_extent; + + /** + * Table for interpolating openlem surface velocities back to ASPECT. + */ + std::array table_intervals; + + /** + * Whether or not to use the ghost nodes. + */ + bool use_ghost_nodes; + + /** + * Check whether openlem needs to be restarted. This is used as + * a mutable bool because we determine whether the model is being resumed in + * initialize(), and then after reinitializing openlem we change it to false + * so it does not initialize openlem again in future timesteps. + * TODO: There is probably a better way to do this, and restarts should be rolled into + * the general ASPECT restart. + */ + mutable bool restart; + + /** + * How many levels openlem should be refined above the maximum ASPECT surface resolution. + */ + unsigned int additional_refinement_levels; + + /** + * Maximum expected refinement level at ASPECT's surface. + * This and resolution_difference are required to properly transfer node data from + * ASPECT to openlem. + */ + unsigned int maximum_surface_refinement_level; + + /** + * User set openlem Y extent for a 2D ASPECT model. + */ + double openlem_y_extent_2d; + + /** + * 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; + + /** + * Suggestion for the number of openlem steps to run for every ASPECT timestep, + * where the openlem timestep is determined by ASPECT_timestep_length divided by + * this parameter. + */ + unsigned int openlem_steps_per_aspect_step; + + /** + * Maximum timestep allowed for openlem, if the suggested timestep exceeds this + * limit it is repeatedly divided by 2 until the final timestep is smaller than this parameter. + */ + double maximum_openlem_timestep; + + /** + * Difference in refinement levels expected at the ASPECT surface, + * where this would be set to 2 if 3 refinement levels are set at the surface. + * This and surface_resolution are required to properly transfer node data from + * ASPECT to FastScape. + * + * TODO: Should this be kept this way, or make it so the input is the expected levels + * of refinement at the surface, and we can subtract one within the code? Also, + * it would be good to find a way to check these are correct, because they are a + * common source of errors. + */ + unsigned int surface_refinement_difference; + + + /** + * Node tolerance for how close a ASPECT node must be to the FastScape node + * for the value to be transferred. This is only necessary if use_v is set to 0 + * and the free surface is used to advect the surface with a normal projection, or + * if there is a surface refinement level difference leading to excess interpolation + * points in areas of high ASPECT resolution. + */ + double node_tolerance; + + /** + * The velocity interpolation data + */ + Table velocity_table; + }; + } +} + + +#endif diff --git a/source/mesh_deformation/openlem.cc b/source/mesh_deformation/openlem.cc new file mode 100644 index 00000000000..6e3f03acbbc --- /dev/null +++ b/source/mesh_deformation/openlem.cc @@ -0,0 +1,1924 @@ +/* + Copyright (C) 2018 - 2023 by the authors of the ASPECT code. + + This file is part of ASPECT. + + ASPECT is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2, or (at your option) + any later version. + + ASPECT is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ASPECT; see the file LICENSE. If not see + . + */ + + + +#include "aspect/global.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace aspect +{ + namespace MeshDeformation + { + template + OpenLEM::OpenLEM() + : + connector(&grid_new) + {} + + + template + void + OpenLEM::initialize () + { + //if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + CitationInfo::add("openlem"); + + AssertThrow(Plugins::plugin_type_matches>(this->get_geometry_model()), + ExcMessage("OpenLEM can only be run with a box geometry model.")); + + const GeometryModel::Box *geometry + = dynamic_cast*> (&this->get_geometry_model()); + + // Find the id associated with the top boundary and boundaries that call mesh deformation. + const types::boundary_id top_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"); + const std::set mesh_deformation_boundary_ids + = this->get_mesh_deformation_handler().get_active_mesh_deformation_boundary_indicators(); + + // Get the deformation type names called for each boundary. + std::map> mesh_deformation_boundary_indicators_map + = this->get_mesh_deformation_handler().get_active_mesh_deformation_names(); + + // Loop over each mesh deformation boundary, and make sure openlem is only called on the surface. + for (const types::boundary_id id : mesh_deformation_boundary_ids) + { + const std::vector &names = mesh_deformation_boundary_indicators_map[id]; + for (const auto &name : names) + { + if (name == "openlem") + AssertThrow(id == top_boundary, + ExcMessage("OpenLEM can only be called on the surface boundary.")); + } + } + + // Several compositional fields are commonly used in conjunction with the openlem plugin, i.e. + // "sediment_age" to track the age of the sediment deposited and "deposition_depth" to track the depth + // with respect to the unperturbed surface of the model domain. Their values are controlled by setting + // boundary conditions on the top boundary that is deformed by openlem. While it is useful to track these + // fields, they are not needed for any function in the openlem plugin. If they exist however, we need + // to make sure that these fields do not have the type "chemical composition" and are therefore not taken + // into account when computing material properties. + const std::vector chemical_field_names = this->introspection().chemical_composition_field_names(); + if (this->introspection().compositional_name_exists("sediment_age")) + { + const std::vector::const_iterator + it = std::find(chemical_field_names.begin(), chemical_field_names.end(), "sediment_age"); + AssertThrow (it == chemical_field_names.end(), + ExcMessage("There is a field sediment_age that is of type chemical composition. " + "Please change it to type generic so that it does not affect material properties.")); + } + if (this->introspection().compositional_name_exists("deposition_depth")) + { + const std::vector::const_iterator + it = std::find(chemical_field_names.begin(), chemical_field_names.end(), "deposition_depth"); + AssertThrow (it == chemical_field_names.end(), + ExcMessage("There is a field deposition_depth that is of type chemical composition. " + "Please change it to type generic so that it does not affect material properties.")); + } + + // Initialize parameters for restarting openlem + restart = this->get_parameters().resume_computation; + + // Since we don't open these until we're on one process, we need to check if the + // restart files exist beforehand. + if (restart) + { + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + AssertThrow(Utilities::fexists(this->get_output_directory() + "openlem_elevation_restart.txt"), + ExcMessage("Cannot open topography file to restart openLEM.")); + AssertThrow(Utilities::fexists(this->get_output_directory() + "openlem_basement_restart.txt"), + ExcMessage("Cannot open basement file to restart openLEM.")); + AssertThrow(Utilities::fexists(this->get_output_directory() + "openLEM_silt_fraction_restart.txt"), + ExcMessage("Cannot open silt fraction file to restart openLEM.")); + } + } + // The first entry represents the minimum coordinates of the model domain, the second the model extent. + for (unsigned int d=0; dget_origin()[d]; + grid_extent[d].second = geometry->get_extents()[d]; + } + + // Get the x and y repetitions used in the parameter file so + // the openlem cell size can be properly set. + const std::array repetitions = geometry->get_repetitions(); + + // Set number of x points, which is generally 1+(openlem refinement level)^2. + // The openlem refinement level is a combination of the maximum ASPECT refinement level + // at the surface and any additional refinement we want in openlem. If + // repetitions are specified we need to adjust the number of points to match what ASPECT has, + // which can be determined by multiplying the points by the repetitions before adding 1. + // Finally, if ghost nodes are used we add two additional points on each side. + const unsigned int multiply_openlem_extent_factor = 2; + const unsigned int ghost_nodes = 2*use_ghost_nodes; + const unsigned int openlem_refinement_level = maximum_surface_refinement_level + additional_refinement_levels; + const unsigned int aspect_refinement_level = maximum_surface_refinement_level; + const unsigned int openlem_nodes = Utilities::pow(2,openlem_refinement_level); + const unsigned int aspect_nodes = Utilities::pow(2,aspect_refinement_level); + openlem_nx = openlem_nodes * repetitions[0] + ghost_nodes + 1; + aspect_nx = aspect_nodes * repetitions[0] + ghost_nodes + 1; + + // Size of openlem cell. + openlem_dx = (grid_extent[0].second)/(openlem_nodes * repetitions[0]); + aspect_dx = (grid_extent[0].second)/(aspect_nodes * repetitions[0]); + + // openlem X extent, which is generally ASPECT's extent unless the ghost nodes are used, + // in which case 2 cells are added on either side. + openlem_x_extent = (grid_extent[0].second* multiply_openlem_extent_factor) + openlem_dx * ghost_nodes; + aspect_x_extent = (grid_extent[0].second) + aspect_dx * ghost_nodes; + + // Sub intervals are 3 less than points, if including the ghost nodes. Otherwise 1 less. + //table_intervals[0] = openlem_nodes * repetitions[0]; + //table_intervals[dim-1] = 1; + + if (dim == 2) + { + openlem_dy = openlem_dx; + openlem_y_extent = std::round(openlem_y_extent_2d/openlem_dy)*openlem_dy + openlem_dy * ghost_nodes; + openlem_ny = 1+(openlem_y_extent/openlem_dy); + aspect_dy = aspect_dx; + aspect_y_extent = std::round(openlem_y_extent_2d/aspect_dy)*aspect_dy + aspect_dy * ghost_nodes; + aspect_ny = 1+aspect_y_extent/aspect_dy; + } + else + { + openlem_ny = openlem_nodes * repetitions[1] + ghost_nodes + 1; + openlem_dy = (grid_extent[1].second)/(openlem_nodes * repetitions[1]); + //table_intervals[1] = openlem_nodes * repetitions[1]; + openlem_y_extent = (grid_extent[1].second * multiply_openlem_extent_factor) + openlem_dy *ghost_nodes; + aspect_ny = aspect_nodes * repetitions[1] + ghost_nodes + 1; + aspect_dy = (grid_extent[1].second)/(aspect_nodes * repetitions[1]); + aspect_y_extent = (grid_extent[1].second) + aspect_dy * ghost_nodes; + } + + openlem_x_extent *= multiply_openlem_extent_factor; + openlem_y_extent *= multiply_openlem_extent_factor; + openlem_nx *= multiply_openlem_extent_factor; + openlem_ny *= multiply_openlem_extent_factor; + + // Create a folder for the openlem visualization files. + Utilities::create_directory (this->get_output_directory() + "openLEM/", + this->get_mpi_communicator(), + false); + + last_output_time = 0; + // preparte the openLEM variables + grid_old = openlem::OceanGrid(openlem_nx,openlem_ny); + grid_new = openlem::OceanGrid(openlem_nx,openlem_ny); + std::cout << "nx = " << openlem_nx << ", ny = " << openlem_ny << ", grid_extent.first = " << grid_extent[0].first << ":" << grid_extent[1].first<< ", grid_extent.second = " << grid_extent[0].second << ":" << grid_extent[1].second << ", openlem extent = " << openlem_x_extent << ":" << openlem_y_extent << std::endl; + + connector = openlem::Connector(&grid_new,openlem_dy,grid_extent[0].first-(grid_extent[0].second/multiply_openlem_extent_factor), grid_extent[1].first-(grid_extent[1].second/multiply_openlem_extent_factor)); + // todo: fill the grid new + + // Define all nodes with non-positive elevations (here, the ocean around the + // island as boundary nodes + std::mt19937 random_number_generator(1);//(openlem_seed); + std::uniform_real_distribution random_distribution(-1.0,1.0);//(-noise_elevation,noise_elevation); + for ( int i = 0; i < grid_new.m; ++i ) + for ( int j = 0; j < grid_new.n; ++j ) + { + //grid_new[i][j].b = grid_new[i][j].h <= 0 || i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n-1; + grid_new[i][j].b = i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n-1 + || i == 1 || j == 1 || i == grid_new.m-2 || j == grid_new.n-2 + ; + //grid_new[i][j].h = (grid_new[i][j].b = + // i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n-1) ? 1 : this->get_initial_topography_model().value(Point(connector.x[i][j],connector.y[i][j])); //|| + //i == 1 || j == 1 || i == grid_new.m-2 || j == grid_new.n-2 + //; + grid_new[i][j].h = i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n-1 || i == 1 || j == 1 || i == grid_new.m-2 || j == grid_new.n-2 ? 1 : this->get_initial_topography_model().value(Point(connector.x[i][j],connector.y[i][j]))+random_distribution(random_number_generator); + //grid_new[i][j].b += (grid_new[i][j].h <= 0)*2; + } + + // Compute initial flow pattern and water level + deepest_point = openlem::Point(0,0); + for ( int i = 0; i < grid_new.m; ++i ) + for ( int j = 0; j < grid_new.n; ++j ) + { + if ( grid_new[i][j].h < grid_new[deepest_point].h ) + deepest_point = openlem::Point(i,j); + } + grid_new[deepest_point].b = 1; + grid_new.odiff = openlem_ocean_diffusivity; + grid_new.fillLakes(); + grid_new.computeFluxes(); + //grid_old = grid_new; + std::cout << "opelem_dy = " << openlem_dy << ", dox = " << grid_extent[0].first << ", doy" << grid_extent[1].first << ", deepest point = " << deepest_point.i << ":" << deepest_point.j << ", q of deepest point = " << grid_new[deepest_point].q << ", deepest point h = " << grid_new[deepest_point].h << ", b = " << (int)grid_new[deepest_point].b << std::endl; + //connector = openlem::Connector(&grid_new,openlem_dy,grid_extent[0].first, grid_extent[1].first, 0.5); + //connector.interpolation(double &x, double &y) + for (unsigned int x_i = 0; x_i < openlem_nx; ++x_i) + for (unsigned int y_i= 0; y_i < openlem_ny; ++y_i) + { + //std::cout << "a x:y = " << x_i << ":" << y_i << " = " << connector.x[x_i][y_i] << ":" << connector.y[x_i][y_i] << std::endl; + //connector.interpolation(connector.x[x_i][y_i], connector.y[x_i][y_i]); + } + //connector = openlem::Connector(&grid_new,openlem_dy,grid_extent[0].first, grid_extent[1].first ); + //for (unsigned int x_i = 0; x_i < openlem_nx; ++x_i) + // for (unsigned int y_i= 0; y_i < openlem_ny; ++y_i) + // std::cout << "x:y = " << x_i << ":" << y_i << " = " << connector.x[x_i][y_i] << ":" << connector.y[x_i][y_i] << std::endl; + } + openlem_mesh_aspect_cell = std::vector>(openlem_nx,std::vector(openlem_ny)); + aspect_mesh_aspect_cell_id = std::vector>(aspect_nx,std::vector(aspect_ny)); + aspect_mesh_aspect_cell_face_number = std::vector>(aspect_nx,std::vector(aspect_ny,std::numeric_limits::signaling_NaN())); + //for (unsigned int xi = 0; xi < aspect_nx; ++xi) + // { + // for (unsigned int yi = 0; yi < aspect_ny; ++yi) + // { + // const double x_coord = xi*aspect_dx + 0.5*aspect_dx - grid_extent[0].first; + // const double y_coord = yi*aspect_dy + 0.5*aspect_dy - grid_extent[1].first; + // aspect_mesh_aspect_cell[xi][yi] = 1; + // } + // } + + + } + + + template + void + OpenLEM::update () + { + if (this->get_timestep_number() == 0) + { + aspect_mesh_aspect_cell_id = std::vector>(aspect_nx,std::vector(aspect_ny)); + aspect_mesh_aspect_cell_face_number = std::vector>(aspect_nx,std::vector(aspect_ny,std::numeric_limits::signaling_NaN())); + + const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"); + // Get a quadrature rule that exists only on the corners, and increase the refinement if specified. + const QIterated face_corners (QTrapezoid<1>(), + 1);//additional_refinement_levels+surface_refinement_difference)); + + FEFaceValues fe_face_values (this->get_mapping(), + this->get_fe(), + face_corners, + update_values | + update_quadrature_points); + // first check if the stored cell is still the cell we need. + //std::cout << "face_corners.size() = " << face_corners.size() << std::endl; + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + if (cell->is_locally_owned() && cell->at_boundary()) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) + if (cell->face(face_no)->at_boundary()) + { + if ( cell->face(face_no)->boundary_id() != relevant_boundary) + continue; + + fe_face_values.reinit(cell, face_no); + const Point vertex = fe_face_values.get_cell()->center(); + const double x_coord = vertex[0] - grid_extent[0].first; + const double y_coord = vertex[1] - grid_extent[1].first; + const size_t ixa = std::floor(x_coord/aspect_dx); + const size_t iya = std::floor(y_coord/aspect_dy); + //std::cout << "x:y = " << ixa << ":" << iya << ", vertex = " << vertex << std::endl; + aspect_mesh_aspect_cell_id[ixa][iya] = cell->id(); + aspect_mesh_aspect_cell_face_number[ixa][iya] = face_no; + } + + for (unsigned int i = 0; i < aspect_nx-1; ++i) + { + aspect_mesh_aspect_cell_id[i][aspect_ny-1] = aspect_mesh_aspect_cell_id[i][aspect_ny-2]; + aspect_mesh_aspect_cell_face_number[i][aspect_ny-1] = aspect_mesh_aspect_cell_face_number[i][aspect_ny-2]; + } + for (unsigned int i = 0; i < aspect_ny; ++i) + { + aspect_mesh_aspect_cell_id[aspect_nx-1][i] = aspect_mesh_aspect_cell_id[aspect_nx-2][i]; + aspect_mesh_aspect_cell_face_number[aspect_nx-1][i] = aspect_mesh_aspect_cell_face_number[aspect_nx-2][i]; + } + for (unsigned int xia = 0; xia < aspect_nx; ++xia) + { + for (unsigned int yia = 0; yia < aspect_nx; ++yia) + { + if (aspect_mesh_aspect_cell_id[xia][yia].get_coarse_cell_id() != numbers::invalid_coarse_cell_id) + { + const auto &cell = this->get_triangulation().create_cell_iterator(aspect_mesh_aspect_cell_id[xia][yia])->as_dof_handler_iterator(this->get_dof_handler()); + + fe_face_values.reinit(cell, aspect_mesh_aspect_cell_face_number[xia][yia]); + const Point vertex = fe_face_values.get_cell()->center(); + //std::cout << "x:y = " << xia << ":" << yia << ", fn: = " << aspect_mesh_aspect_cell_face_number[xia][yia] << ", cellid = " << aspect_mesh_aspect_cell_id[xia][yia].get_coarse_cell_id() << ", 2d cell center = " << vertex << std::endl; + } + } + } + } + std::string dirname = this->get_output_directory(); + const char *dirname_char=dirname.c_str(); + const unsigned int dirname_length = dirname.length(); + + TimerOutput::Scope timer_section(this->get_computing_timer(), "openlem plugin"); + + std::cout << "Flag 1: before (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + if (this->get_timestep_number() == 0) + { + for (unsigned int x_i=0; x_ih = this->get_initial_topography_model().value(Point(connector.x[x_i][y_i],connector.y[x_i][y_i])); + } + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + connector.write_vtk(grid_extent[dim-1].second, 0, 0., dirname); + return; + } + grid_old = grid_new; + //connector = openlem::Connector(&grid_new,openlem_dy); + std::cout << "Flag 2 before (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + + const unsigned int current_timestep = this->get_timestep_number (); + const double aspect_timestep_in_years = this->get_timestep() / year_in_seconds; + + // Find a openlem timestep that is below our maximum timestep. + unsigned int openlem_iterations = openlem_steps_per_aspect_step; + double openlem_timestep_in_years = aspect_timestep_in_years/openlem_iterations; + while (openlem_timestep_in_years>maximum_openlem_timestep) + { + openlem_iterations *= 2; + openlem_timestep_in_years *= 0.5; + } + + // Vector to hold the velocities that represent the change to the surface. + //const unsigned int openlem_array_size = openlem_nx*openlem_ny; + std::vector>> mesh_locations(openlem_nx,std::vector>(openlem_ny)); + // fill mesh locations + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + for (unsigned int x_i = 0; x_i < openlem_nx; ++x_i) + { + for (unsigned int y_i = 0; y_i < openlem_nx; ++y_i) + { + // The z will be replaced later on + if (dim == 3) + { + mesh_locations[x_i][y_i] = Point(connector.x[x_i][y_i],connector.y[x_i][y_i],grid_extent[dim-1].second + grid_new[x_i][y_i].h); + //if (x_i == 10) + // std::cout << "A: rank = " <get_mpi_communicator()) << ", x:y = " << x_i << ":" << y_i << ", mesh_locations[x_i][y_i] = " << mesh_locations[x_i][y_i] << ", connector x:y = " << connector.x[x_i][y_i] << ":" << connector.x[x_i][y_i] << std::endl; + // NOTE: at timestep 0 + //if (this->get_timestep_number() == 0) + // { + // // TODO: probably should move this if statement out of the double loop + // //std::cout << "it: " << this->get_initial_topography_model().value(Point(connector.x[x_i][y_i],connector.y[x_i][y_i])) << ", x = " << connector.x[x_i][y_i] << ", y = " << connector.y[x_i][y_i] << std::endl; + // mesh_locations[x_i][y_i] = Point(connector.x[x_i][y_i],connector.y[x_i][y_i],grid_extent[dim-1].second);// + this->get_initial_topography_model().value(Point(connector.x[x_i][y_i],connector.y[x_i][y_i]))); + // } + //else + // { + // mesh_locations[x_i][y_i] = Point(connector.x[x_i][y_i],connector.y[x_i][y_i],grid_extent[dim-1].second);// + grid_new[x_i][y_i].h); + // } + } + } + } + } + + + std::cout << "before broadcast mesh_locations" << std::endl; + mesh_locations = Utilities::MPI::broadcast(this->get_mpi_communicator(),mesh_locations,0); + std::cout << "after broadcast mesh_locations" << std::endl; + + Tensor<1,dim> invalid_velocity; + for (unsigned int dim_i = 0; dim_i < dim-1; ++dim_i) + { + invalid_velocity[dim_i] = std::numeric_limits::quiet_NaN(); + } + invalid_velocity[dim-1] = 0; + std::vector>> openlem_mesh_velocities(openlem_nx,std::vector>(openlem_ny, invalid_velocity)); + //std::vector>> mesh_velocity_locations(openlem_nx,std::vector>(openlem_ny)); + //std::vector> mesh_velocity_distances(openlem_nx,std::vector(openlem_ny)); + aspect_mesh_z = std::vector>(aspect_nx,std::vector(aspect_ny,std::numeric_limits::signaling_NaN())); + + + std::vector>>> mesh_velocities_vector; + + //const UpdateFlags update_flags = update_values | update_gradients; + //std::unique_ptr> evaluator = construct_solution_evaluator(*this,update_flags); + //std::vector evaluation_flags (this->introspection().n_components, EvaluationFlags::nothing); + + //for (unsigned int i=0; iintrospection().n_components; ++i) + // { + // evaluation_flags[i] |= EvaluationFlags::values; + // evaluation_flags[i] |= EvaluationFlags::gradients; + // } + + + + const UpdateFlags update_flags = update_values;// | update_gradients; + std::unique_ptr> evaluator = construct_solution_evaluator(*this,update_flags); + + // Todo: assert that vectors of cells, positions and reference positino are the same size + + // FEPointEvaluation uses different evaluation flags than the common UpdateFlags. + // Translate between the two. + std::vector evaluation_flags (this->introspection().n_components, EvaluationFlags::nothing); + + for (unsigned int i=0; iintrospection().n_components; ++i) + { + evaluation_flags[i] |= EvaluationFlags::values; + //evaluation_flags[i] |= EvaluationFlags::gradients; + } + + std::vector> velocity_field; + std::vector> solution(this->get_fe().dofs_per_cell); + + //std::vector>> gradients(this->get_fe().dofs_per_cell); + solution.resize(1,std::vector(evaluator->n_components(), numbers::signaling_nan())); + //gradients.resize(1,std::vector>(evaluator->n_components(), numbers::signaling_nan>())); + + const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"); + // Get a quadrature rule that exists only on the corners, and increase the refinement if specified. + const QIterated face_corners (QTrapezoid<1>(), + 1);//additional_refinement_levels+surface_refinement_difference)); + + FEFaceValues fe_face_values (this->get_mapping(), + this->get_fe(), + face_corners, + update_values | + update_quadrature_points); + + const CellId unitialized_cell_id = CellId(); + //Tensor<1,dim> invalid_velocity; + //for (unsigned int dim_i = 0; dim_i < dim; ++dim_i) + // { + // velocity[dim_i] = std::numeric_limits::signaling_NaN(); + // } + for (unsigned int xi = 0; xi < openlem_nx; ++xi) + { + for (unsigned int yi = 0; yi < openlem_ny; ++yi) + { + // get corresponding aspect index + Point<2> openlem_point_2d(connector.x[xi][yi],connector.y[xi][yi]); + + double x_coord = connector.x[xi][yi] - grid_extent[0].first; + double y_coord = connector.y[xi][yi] - grid_extent[1].first; + if ( + x_coord >= 0 && x_coord <= grid_extent[0].second && + y_coord >= 0 && y_coord <= grid_extent[1].second + ) + { + const double x_coord = openlem_point_2d[0] - grid_extent[0].first; + const double y_coord = openlem_point_2d[1] - grid_extent[1].first; + const size_t xia = std::floor(x_coord/aspect_dx); + const size_t yia = std::floor(y_coord/aspect_dy); + //if (xi < openlem_nx-2 && yi > 30 && yi < 50) + // std::cout << "F: openlem_point_2d = " << openlem_point_2d << ", x:ycoord = " << x_coord << ":" << y_coord << ", x:ycoord/dx:dy = " << x_coord/aspect_dx << ":" << y_coord/aspect_dy << ", xia:yia = " << xia << ":" << yia << ", aspect dx:dy = " << aspect_dx << ":" << aspect_dy << std::endl; + + if (aspect_mesh_aspect_cell_id[xia][yia].get_coarse_cell_id() != numbers::invalid_coarse_cell_id) + { + if (!std::isnan(aspect_mesh_aspect_cell_face_number[xia][yia])) + { + const auto &cell = this->get_triangulation().create_cell_iterator(aspect_mesh_aspect_cell_id[xia][yia])->as_dof_handler_iterator(this->get_dof_handler()); + if (cell->state() == IteratorState::valid && cell->is_locally_owned() ) + { + bool found_a_value = false; + std::vector> unit_point(1); + double closest_distance_cell = std::numeric_limits::infinity(); + fe_face_values.reinit(cell, aspect_mesh_aspect_cell_face_number[xia][yia]); + for (unsigned int corner = 0; corner < face_corners.size(); ++corner) + { + const Point vertex = fe_face_values.quadrature_point(corner); + + const double x_coord = vertex[0] - grid_extent[0].first; + const double y_coord = vertex[1] - grid_extent[1].first; + const size_t ixa = std::round(x_coord/aspect_dx); + const size_t iya = std::round(y_coord/aspect_dy); + aspect_mesh_z[ixa][iya] = vertex[dim-1] - grid_extent[dim-1].second; + //mesh_locations[ixa][iya][dim-1] = vertex[dim-1]; + } + { + + this->get_mapping().transform_points_real_to_unit_cell(cell, {mesh_locations[xi][yi]},unit_point); + const double distance_to_unit_cell = GeometryInfo::distance_to_unit_cell(unit_point[0]); + + //if (xi < openlem_nx-2 && yi > 30 && yi < 50) + // std::cout << "U1: velocity " << xi << ":" << yi << ", aspect x:y = " << xia << ":" < cp2d(closest_point_real[0],closest_point_real[1]); + + auto &velocity_evaluator = evaluator->get_velocity_or_fluid_velocity_evaluator(false); + small_vector solution_values(this->get_fe().dofs_per_cell); + cell->get_dof_values(this->get_current_linearization_point(),//input_solution, + solution_values.begin(), + solution_values.end()); + + + solution[0] = std::vector(evaluator->n_components(), numbers::signaling_nan()); + + + evaluator->reinit(cell, closest_point); + //evaluator->evaluate({solution_values.data(),solution_values.size()},evaluation_flags); + //evaluator->get_solution(0, {&solution[0][0],solution[0].size()}, evaluation_flags); + + velocity_evaluator.evaluate({solution_values.data(),solution_values.size()}, + EvaluationFlags::values); + + auto &mapping_info = evaluator->get_mapping_info(); + mapping_info.reinit(cell, {closest_point.data(),closest_point.size()}); + + + //Point<3> velocity_analytical(-openlem_point_2d[1],openlem_point_2d[0],0.0); + Tensor<1,dim> velocity = velocity_evaluator.get_value(0)*year_in_seconds; + //if (xi < openlem_nx-2 && yi > 30 && yi < 50) + // std::cout << "X: velocity " << xi << ":" << yi << " = " << velocity << std::endl; + + if (!isnan(velocity[0])) + { + openlem_mesh_velocities[xi][yi] = velocity; + //if (xi < openlem_nx-2 && yi > 30 && yi < 50) + // std::cout << "Y: velocity " << xi << ":" << yi << " = " << velocity << ", mesh_velocities[xi][yi] = " << openlem_mesh_velocities[xi][yi] << std::endl; + } + } + } + } + } + + } + } + //else + // { + // mesh_velocities[xi][yi] = velocity; + // } + } + } + // first check if the stored cell is still the cell we need. + //std::cout << "face_corners.size() = " << face_corners.size() << std::endl; + /* + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + if (cell->is_locally_owned() && cell->at_boundary()) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) + if (cell->face(face_no)->at_boundary()) + { + if ( cell->face(face_no)->boundary_id() != relevant_boundary) + continue; + + //std::vector> vel(face_corners.size()); + fe_face_values.reinit(cell, face_no); + //fe_face_values[this->introspection().extractors.velocities].get_function_values(this->get_solution(), vel); + + for (unsigned int corner = 0; corner < face_corners.size(); ++corner) + { + const Point vertex = fe_face_values.quadrature_point(corner); + + const double x_coord = vertex[0] - grid_extent[0].first; + const double y_coord = vertex[1] - grid_extent[1].first; + const size_t ixa = std::round(x_coord/aspect_dx); + const size_t iya = std::round(y_coord/aspect_dy); + aspect_mesh_z[ixa][iya] = vertex[dim-1] - grid_extent[dim-1].second; + + } + } + */ + //unsigned int total_found_velocities = 0; + //const CellId unitialized_cell_id = CellId(); + /* + for (unsigned int x_i = 0; x_i < openlem_nx; ++x_i) + { + for (unsigned int y_i = 0; y_i < openlem_ny; ++y_i) + { + bool found_a_value = false; + std::vector> unit_point(1); + double closest_distance_cell = std::numeric_limits::infinity(); + //std::cout << "cellid = " << openlem_mesh_aspect_cell[x_i][y_i].get_coarse_cell_id() << "..." << std::endl; + //std::cout << "try cell id at " << x_i << ":" << y_i << ", coarse cell id = " << openlem_mesh_aspect_cell[x_i][y_i].get_coarse_cell_id()<< std::endl; + if (openlem_mesh_aspect_cell[x_i][y_i].get_coarse_cell_id() != numbers::invalid_coarse_cell_id) + { + //std::cout << "found a valid cell id at " << x_i << ":" << y_i << std::endl; + const auto &cell = this->get_triangulation().create_cell_iterator(openlem_mesh_aspect_cell[x_i][y_i])->as_dof_handler_iterator(this->get_dof_handler()); + if ( cell->state() == IteratorState::valid && cell->is_locally_owned() && cell->at_boundary()) + { + //std::cout << "found a valid cell at " << x_i << ":" << y_i << std::endl; + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) + { + if (cell->face(face_no)->at_boundary()) + { + if ( cell->face(face_no)->boundary_id() != relevant_boundary) + continue; + + this->get_mapping().transform_points_real_to_unit_cell(cell, {mesh_locations[x_i][y_i]},unit_point); + const double distance_to_unit_cell = GeometryInfo::distance_to_unit_cell(unit_point[0]); + + if (distance_to_unit_cell < closest_distance_cell) + { + std::vector> closest_point = {cell->reference_cell().closest_point(unit_point[0])}; + Point closest_point_real = this->get_mapping().transform_unit_to_real_cell(cell,closest_point[0]); + + Point<2> cp2d(closest_point_real[0],closest_point_real[1]); + Point<2> openlem_point_2d(connector.x[x_i][y_i],connector.y[x_i][y_i]); + + + if (cell->state() != IteratorState::valid) + { + continue; + } + + auto &velocity_evaluator = evaluator->get_velocity_or_fluid_velocity_evaluator(false); + small_vector solution_values(this->get_fe().dofs_per_cell); + cell->get_dof_values(this->get_current_linearization_point(),//input_solution, + solution_values.begin(), + solution_values.end()); + + + solution[0] = std::vector(evaluator->n_components(), numbers::signaling_nan()); + + + evaluator->reinit(cell, closest_point); + //evaluator->evaluate({solution_values.data(),solution_values.size()},evaluation_flags); + //evaluator->get_solution(0, {&solution[0][0],solution[0].size()}, evaluation_flags); + + velocity_evaluator.evaluate({solution_values.data(),solution_values.size()}, + EvaluationFlags::values); + + auto &mapping_info = evaluator->get_mapping_info(); + mapping_info.reinit(cell, {closest_point.data(),closest_point.size()}); + + + Point<3> velocity_3d; + Point<3> velocity_original; + Point<3> velocity_analytical(-openlem_point_2d[1],openlem_point_2d[0],0.0); + Tensor<1,dim> velocity = velocity_evaluator.get_value(0)*year_in_seconds; + + //if (y_i == 10) + // std::cout << "rank = " <get_mpi_communicator()) << ", x:y = " << x_i << ":" << y_i << ", closest_point = " << closest_point[0] << ", closest_point_real = " << closest_point_real<< ", mesh_locations = " << mesh_locations[x_i][y_i] << ", velo = " << velocity << ", velo ana = " << velocity_analytical << ", distance_to_unit_cell = " << distance_to_unit_cell<< ", vel raw = " << velocity_evaluator.get_value(0) << std::endl; + for (unsigned int i = 0; i < dim; ++i) + { + velocity_original[i] = velocity[i]; + //velocity[i] = velocity_analytical[i];// * year_in_seconds;//solution_values[this->introspection().component_indices.velocities[i]] * year_in_seconds; + velocity_3d[i] = velocity[i];//solution_values[this->introspection().component_indices.velocities[i]] * year_in_seconds; + } + + // TODO: maybe don't need the closest point real in the tuple? + //std::tuple, Tensor<1,dim>> tuple = std::make_tuple(x_i, y_i, closest_point_real, velocity); + if (!isnan(velocity[0])) + { + found_a_value = true; + closest_distance_cell = distance_to_unit_cell; // TODO: probably just skip out if distance is exactly zero. + //total_found_velocities++; + mesh_velocities_vector.emplace_back(std::make_tuple(x_i, y_i, distance_to_unit_cell, closest_point_real, velocity)); + openlem_mesh_aspect_cell[x_i][y_i] = cell->id(); + //std::cout << "set cell id at " << x_i << ":" << y_i << ", coarse cell id = " << cell->id().get_coarse_cell_id() << ", " << openlem_mesh_aspect_cell[x_i][y_i] << std::endl; + break; + } + } + } + } + } + } + if (closest_distance_cell != 0.0) + { + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + { + if (closest_distance_cell == 0.0) + { + // if the point is inside the cell it is defined as exactly zero + break; + } + if (cell->is_locally_owned() && cell->at_boundary()) + { + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) + { + if (cell->face(face_no)->at_boundary()) + { + if ( cell->face(face_no)->boundary_id() != relevant_boundary) + continue; + + this->get_mapping().transform_points_real_to_unit_cell(cell, {mesh_locations[x_i][y_i]},unit_point); + const double distance_to_unit_cell = GeometryInfo::distance_to_unit_cell(unit_point[0]); + + + if (distance_to_unit_cell < closest_distance_cell) + { + std::vector> closest_point = {cell->reference_cell().closest_point(unit_point[0])}; + Point closest_point_real = this->get_mapping().transform_unit_to_real_cell(cell,closest_point[0]); + + Point<2> cp2d(closest_point_real[0],closest_point_real[1]); + Point<2> openlem_point_2d(connector.x[x_i][y_i],connector.y[x_i][y_i]); + + + if (cell->state() != IteratorState::valid) + { + continue; + } + + auto &velocity_evaluator = evaluator->get_velocity_or_fluid_velocity_evaluator(false); + small_vector solution_values(this->get_fe().dofs_per_cell); + cell->get_dof_values(this->get_current_linearization_point(),//input_solution, + solution_values.begin(), + solution_values.end()); + + + solution[0] = std::vector(evaluator->n_components(), numbers::signaling_nan()); + + + evaluator->reinit(cell, closest_point); + //evaluator->evaluate({solution_values.data(),solution_values.size()},evaluation_flags); + //evaluator->get_solution(0, {&solution[0][0],solution[0].size()}, evaluation_flags); + + velocity_evaluator.evaluate({solution_values.data(),solution_values.size()}, + EvaluationFlags::values); + + auto &mapping_info = evaluator->get_mapping_info(); + mapping_info.reinit(cell, {closest_point.data(),closest_point.size()}); + + + Point<3> velocity_3d; + Point<3> velocity_original; + Point<3> velocity_analytical(-openlem_point_2d[1],openlem_point_2d[0],0.0); + Tensor<1,dim> velocity = velocity_evaluator.get_value(0)*year_in_seconds; + + //if (y_i == 10) + // std::cout << "rank = " <get_mpi_communicator()) << ", x:y = " << x_i << ":" << y_i << ", closest_point = " << closest_point[0] << ", closest_point_real = " << closest_point_real<< ", mesh_locations = " << mesh_locations[x_i][y_i] << ", velo = " << velocity << ", velo ana = " << velocity_analytical << ", distance_to_unit_cell = " << distance_to_unit_cell<< ", vel raw = " << velocity_evaluator.get_value(0) << std::endl; + for (unsigned int i = 0; i < dim; ++i) + { + velocity_original[i] = velocity[i]; + //velocity[i] = velocity_analytical[i];// * year_in_seconds;//solution_values[this->introspection().component_indices.velocities[i]] * year_in_seconds; + velocity_3d[i] = velocity[i];//solution_values[this->introspection().component_indices.velocities[i]] * year_in_seconds; + } + + // TODO: maybe don't need the closest point real in the tuple? + //std::tuple, Tensor<1,dim>> tuple = std::make_tuple(x_i, y_i, closest_point_real, velocity); + if (!isnan(velocity[0])) + { + found_a_value = true; + closest_distance_cell = distance_to_unit_cell; // TODO: probably just skip out if distance is exactly zero. + //total_found_velocities++; + mesh_velocities_vector.emplace_back(std::make_tuple(x_i, y_i, distance_to_unit_cell, closest_point_real, velocity)); + openlem_mesh_aspect_cell[x_i][y_i] = cell->id(); + //std::cout << "set cell id at " << x_i << ":" << y_i << ", coarse cell id = " << cell->id().get_coarse_cell_id() << ", " << openlem_mesh_aspect_cell[x_i][y_i] << std::endl; + break; + } + } + } + } + } + } + } + else + { + //std::cout << "Found the right cell in storage! " << std::endl; + } + AssertThrow(found_a_value,ExcMessage("Could not found a value for " + std::to_string(x_i) + ":" + std::to_string(y_i))); + } + } + */ + + // the mesh velocties vector should contain all the entries now, return to process 0 and unpack + std::cout << "before gather mesh_velocites_vector" << std::endl; + std::vector>> aspect_mesh_z_vectors = Utilities::MPI::gather(this->get_mpi_communicator(),aspect_mesh_z,0); + //std::vector, Tensor<1,dim>>>> mesh_velocties_vectors = Utilities::MPI::gather(this->get_mpi_communicator(),mesh_velocities_vector,0); + std::vector>>> mesh_velocties_vectors = Utilities::MPI::gather(this->get_mpi_communicator(),openlem_mesh_velocities,0); + std::cout << "after gather mesh_velocites_vector" << std::endl; + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + //std::cout << " aspect_mesh_z_vectors = " << aspect_mesh_z_vectors.size() << std::endl; + for (auto &aspect_mesh_z_item : aspect_mesh_z_vectors) + { + for (size_t xi = 0; xi < aspect_nx; ++xi) + { + for (size_t yi = 0; yi < aspect_ny; ++yi) + { + + //if (xi == 54/4 && yi > 125/4 && yi < 130/4) + //if (xi == 2 && yi > 0 && yi < 10) + // std::cout << "Flag 100: mesh z " << xi << ":" << yi << " = " << aspect_mesh_z_item[xi][yi] << std::endl; + if (!std::isnan(aspect_mesh_z_item[xi][yi])) + { + aspect_mesh_z[xi][yi] = aspect_mesh_z_item[xi][yi]; + //if (xi == 54/4 && yi > 125/4 && yi < 130/4) + //if (xi == 2 && yi > 0 && yi < 10) + // std::cout << "Flag 101: mesh z " << xi << ":" << yi << " = " << aspect_mesh_z[xi][yi] << std::endl; + } + } + + } + } + //std::cout << "mesh_velocties_vectors = " << mesh_velocties_vectors.size() << std::endl; + for (auto &openlem_mesh_velocities_vector : mesh_velocties_vectors) + { + //std::cout << "B: mesh_velocties_vectors = " << mesh_velocties_vectors.size() << std::endl; + for (size_t xi = 0; xi < openlem_nx; ++xi) + { + for (size_t yi = 0; yi < openlem_ny; ++yi) + { + //if (std::get<0>(position_velocity) == 10) + // std::cout << counter << "a: rank = " <get_mpi_communicator()) << ", mesh velocties x:y = " << std::get<0>(position_velocity) << ":" << std::get<1>(position_velocity) << ", velo = " << std::get<3>(position_velocity) << ", distance = " << std::get<2>(position_velocity)<< ", mesh_velocity_distances = " << mesh_velocity_distances[std::get<0>(position_velocity)][std::get<1>(position_velocity)] << std::endl; + //if (xi == 2 && yi > 0 && yi < 10) + // std::cout << "Flag 102: openlem_mesh_velocities " << xi << ":" << yi << " = " << openlem_mesh_velocities_vector[xi][yi] << std::endl; + if (!std::isnan(openlem_mesh_velocities_vector[xi][yi][0])) + { + // if (std::get<0>(position_velocity) == 10) + // std::cout << counter << "b: rank = " <get_mpi_communicator()) << ", mesh velocties x:y = " << std::get<0>(position_velocity) << ":" << std::get<1>(position_velocity) << ", velo = " << std::get<3>(position_velocity) << ", distance = " << std::get<2>(position_velocity)<< std::endl; + openlem_mesh_velocities[xi][yi] = openlem_mesh_velocities_vector[xi][yi]; + //if (xi == 2 && yi > 0 && yi < 10) + // std::cout << "Flag 103: openlem_mesh_velocities " << xi << ":" << yi << " = " << openlem_mesh_velocities[xi][yi] << std::endl; + //mesh_velocity_locations[std::get<0>(position_velocity)][std::get<1>(position_velocity)] = std::get<3>(position_velocity); + //mesh_velocity_distances[std::get<0>(position_velocity)][std::get<1>(position_velocity)] = std::get<2>(position_velocity); + } + } + } + } + // TODO: check whether the x and y of the closest_point_real and the original connector point are actually the same (or close enought) + + + //std::vector> local_aspect_values = get_aspect_values(); + // Because there is no increase in time during timestep 0, we return and only + // initialize and run openlem from timestep 1 and on. + for (unsigned int x_i=0; x_iget_timestep_number() == 0) + { + node->h = this->get_initial_topography_model().value(Point(connector.x[x_i][y_i],connector.y[x_i][y_i])); + } + + node->l = 0; + node->u = openlem_mesh_velocities[x_i][y_i][dim-1];//*100000;//1;//local_aspect_values[dim+2][i]; + // TODO: 2D + connector.vx[x_i][y_i] = openlem_mesh_velocities[x_i][y_i][0];//*100000; + connector.vy[x_i][y_i] = openlem_mesh_velocities[x_i][y_i][1];//*100000; + connector.vz[x_i][y_i] = openlem_mesh_velocities[x_i][y_i][2];//*100000; + } + //connector.write_vtk(grid_extent[dim-1].second,this->get_timestep_number(), 0.0, dirname , "_preconvert"); + connector.convertVelocities(openlem_minimize_advection); + //connector.write_vtk(grid_extent[dim-1].second,this->get_timestep_number(), 0.0, dirname , "_postconvert"); + //connector.write_vtk(grid_extent[dim-1].second, this->get_timestep_number() , "_postconvert"); + connector.computeUpliftRate(); + //connector.write_vtk(grid_extent[dim-1].second,this->get_timestep_number(), 0.0, dirname , "_postuplift"); + //connector.write_vtk(grid_extent[dim-1].second, this->get_timestep_number() , "_postuplift"); + // Define all nodes with non-positive elevations (here, the ocean around the + // island as boundary nodes + for ( int i = 0; i < grid_new.m; ++i ) + for ( int j = 0; j < grid_new.n; ++j ) + { + //grid_new[i][j].b = grid_new[i][j].h <= 0; + // grid_new[i][j].b = i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n; + //grid_new[i][j].b = 0; + //i == 0 || j == 0 || i == grid_new.m-1 || j == grid_new.n-1 || + //i == 1 || j == 1 || i == grid_new.m-2 || j == grid_new.n-2 + //; + //grid_new[i][j].b += (grid_new[i][j].h <= 0)*2; + // grid_new[i][j].u = grid_new[i][j].h <= 0; + } + openlem::OceanGrid &g = grid_new; + grid_new.kt[0] = openlem_kt;//1e-14; + grid_new.kt[1] = openlem_kt;//1e-14; + grid_new.kd[0] = 1e6*openlem_kd;// 1e-14; + grid_new.kd[1] = openlem_kd;// 1e-14; + + //grid_new.fillLakes(); + // find deepest point + // call set deepest point to 1 + // fill lakes + // in each timestep compute water level() + // call function clear ocean() + // call function mark ocean () + // + + } + if (this->get_timestep_number() == 0) + { + connector.write_vtk(grid_extent[dim-1].second, 0, 0.0, dirname); + return; + } + + + // Run openlem on single process. + //mesh_velocity_z = std::vector>(openlem_nx,std::vector(openlem_ny)); + aspect_mesh_dh = std::vector>(aspect_nx,std::vector(aspect_ny)); + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 0) + { + // Initialize the variables that will be sent to openlem. + // Elevation is initialized at a very high number so that we can later check that all points + // received data from ASPECT, and if not throw an assert. + + //std::cout << "Flag 3 before (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + + //for (unsigned int ix=0; ixh = local_aspect_values[0][i]; + // node->l = 0; + // node->u = ;//1;//local_aspect_values[dim+2][i]; + // } + //fill_openlem_arrays(grid_new, + // local_aspect_values); + //std::cout << "Flag 4 before (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + + if (current_timestep == 1 || restart) + { + this->get_pcout() << " Initializing openlem... " << (1+maximum_surface_refinement_level+additional_refinement_levels) << + " levels, cell size: " << openlem_dx << " m." << std::endl; + + // TODO + //initialize_openlem(elevation, + // basement, + // bedrock_transport_coefficient_array, + // bedrock_river_incision_rate_array, + // silt_fraction); + } + else + { + // If it isn't the first timestep we ignore initialization and instead copy all height values from openlem. + // Generally, we overwrite the topography data from ASPECT as openlem may be at a higher resolution. However, + // if we are not using openlem to advect then we do not want to do this and instead use the ASPECT values. + // if (openlem_advection_uplift) + // openlem_copy_h_(elevation.data()); + } + + // Find the appropriate sediment rain based off the time interval. + const double time_in_years = this->get_time();// / year_in_seconds; + //auto it = std::lower_bound(sediment_rain_times.begin(), sediment_rain_times.end(), time_in_years); + //const unsigned int inds = std::distance(sediment_rain_times.begin(), it); + //const double sediment_rain = sediment_rain_rates[inds]; + + //grid_old = grid_new; + std::cout << "main before (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + execute_openlem(grid_new,openlem_timestep_in_years,openlem_iterations,dirname); + + //connector.write_vtk(grid_extent[dim-1].second,this->get_timestep_number(), 0.0, dirname , "_postexectuteopenlem"); + //connector.updateCoordinateSystem (aspect_timestep_in_years); + std::cout << "main after (old:new)= " << (grid_old.getNode(5,7))->h << ":" << (grid_new.getNode(5,7))->h << std::endl; + // Write a file to store h, b & step for restarting. + // TODO: It would be good to roll this into the general ASPECT checkpointing, + // and when we do this needs to be changed. + //if (((this->get_parameters().checkpoint_time_secs == 0) && + // (this->get_parameters().checkpoint_steps > 0) && + // ((current_timestep + 1) % this->get_parameters().checkpoint_steps == 0)) || + // (this->get_time() == this->get_end_time() && this->get_timestepping_manager().need_checkpoint_on_terminate())) + // { + // //save_restart_files(elevation, + // // basement, + // // silt_fraction); + // } + + // Find out our velocities from the change in height. + // Where mesh_velocity_z is a vector of array size that exists on all processes. + //for (unsigned int ix=0; ixh - node_old->h)/(1.0*this->get_timestep());//aspect_timestep_in_years; + // if ((ix == 18 || ix == 19) && iy == 29) + // std::cout << ix << ":" << iy << " = " << mesh_velocity_z[ix][iy] << ", node_old->h = " << node_old->h << ", node_new->h = " << node_new->h << ", node_new->h - node_old->h = " << node_new->h - node_old->h << std::endl; + // } + + //Utilities::MPI::broadcast(this->get_mpi_communicator(), mesh_velocity_z, 0); + + + // loop over all openlem grid points, find the closest aspect grid point + // add the velocity to the closest grid point and count how many grid points contributed to that point + // then divide the value by the counter. + std::vector> aspect_mesh_counter = std::vector>(aspect_nx,std::vector(aspect_ny)); + //double min_openlem_h = std::numeric_limits::infinity(); + //double max_openlem_h = -std::numeric_limits::infinity(); + for (unsigned int xi=0; xi= 0 && x_coord <= grid_extent[0].second && + y_coord >= 0 && y_coord <= grid_extent[1].second + + //connector.x[ix][iy] >= grid_extent[0].first && connector.x[ix][iy] <= grid_extent[0].first+ grid_extent[0].second && + //connector.y[ix][iy] >= grid_extent[1].first && connector.x[ix][iy] <= grid_extent[1].first+ grid_extent[1].second + ) + { + + //Point<2> openlem_point (connector.x[ix][iy],connector.y[ix][iy]); + //Point<2> aspect_point_round (std::round((connector.x[ix][iy]-grid_extent[0].first)/aspect_dx),std::round((connector.y[ix][iy]-grid_extent[1].first)/aspect_dy)); + // TODO: add assert that points are > 0 + size_t closest_point_ixa = std::round(x_coord/aspect_dx);//aspect_point_round[0];//std::numeric_limits::signaling_NaN(); + size_t closest_point_iya = std::round(y_coord/aspect_dy);//aspect_point_round[1];//std::numeric_limits::signaling_NaN(); + //if (grid_new.getNode(ix,iy)->b == 0) + //min_openlem_h = std::min(min_openlem_h,grid_new.getNode(ix,iy)->h); + //max_openlem_h = std::max(max_openlem_h,grid_new.getNode(ix,iy)->h); + aspect_mesh_dh[closest_point_ixa][closest_point_iya] += grid_new.getNode(xi,yi)->h - aspect_mesh_z[closest_point_ixa][closest_point_iya];// + grid_extent[dim-1].second;//- mesh_velocity_locations[ix][iy][dim-1]+grid_extent[dim-1].second ;// grid_old.getNode(ix,iy)->h;//TODO: check if this is alright or if I need to get the closest node in the old grid sepeartly. //mesh_velocity_z[ix][iy]; + if (xi == 2 && yi > 0 && yi < 10) + std::cout << "Flag 201 x:y " << xi << ":" << yi << ", aspect x:y = " << closest_point_ixa << ":" << closest_point_iya << ", aspect_mesh_dh = " << aspect_mesh_dh[closest_point_ixa][closest_point_iya] << ", grid_new.getNode(xi,yi)->h = " << grid_new.getNode(xi,yi)->h << ", aspect_mesh_z = " << aspect_mesh_z[closest_point_ixa][closest_point_iya] << std::endl; + //if (grid_new.getNode(ix,iy)->h > 10 ) + // std::cout << "rank = " << Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) << ", ix:iy = " << ix << ":" << iy << ", ixa:iya = " << closest_point_ixa << ":" << closest_point_iya << ", h = " << grid_new.getNode(ix,iy)->h << ", aspect_mesh_z =" << aspect_mesh_z[closest_point_ixa][closest_point_iya] << ", result = " << aspect_mesh_dh[closest_point_ixa][closest_point_iya] << std::endl; + aspect_mesh_counter[closest_point_ixa][closest_point_iya] += 1; + } + } + + for (unsigned int ixa=0; ixaget_mpi_communicator()) == 0) + { + //std::cout << "ixa:iya = " << ixa << ":" << iya << ",aspect_mesh_elevation_h = " << aspect_mesh_elevation_h[ixa][iya] << ", counter = " << aspect_mesh_counter[ixa][iya]<< std::endl; + } + //if ( grid_new.getNode(ix,iy)->b == 0) + { + if (aspect_mesh_counter[ixa][iya] == 0) + { + aspect_mesh_dh[ixa][iya] = 0; + } + else + { + AssertThrow(aspect_mesh_counter[ixa][iya] != 0, ExcMessage("Division by zero when averaging aspect mesh velocity z")); + //const double before = aspect_mesh_elevation_h[ixa][iya]; + aspect_mesh_dh[ixa][iya] /= (aspect_mesh_counter[ixa][iya]); + //AssertThrow(aspect_mesh_elevation_h[ixa][iya] <= max_openlem_h && aspect_mesh_elevation_h[ixa][iya] >= min_openlem_h , ExcMessage("wrong!!!! min:max h = " + std::to_string(min_openlem_h) + ":" + std::to_string(max_openlem_h) + ", value:after:counter = " + std::to_string(aspect_mesh_elevation_h[ixa][iya]) + ":" + std::to_string(before) + ":" + std::to_string(aspect_mesh_counter[ixa][iya]))); + } + //if ((ixa == 18*2 || ixa == 19*2) && iya == 29*2) + //std::cout << "ix:iy = " << ixa << ":" << iya << ", aspect_mesh_elevation_h = " << aspect_mesh_elevation_h[ixa][iya] << std::endl; + } + + } + } + else + { + std::cout << "main before (old:new) " << std::endl; + execute_openlem(grid_new,openlem_timestep_in_years,openlem_iterations, dirname); + + std::cout << "main after (old:new) " << std::endl; + //for (unsigned int i=0; iget_mpi_communicator()); + + // Check whether the openlem mesh was filled with data. + // const bool openlem_mesh_filled = Utilities::MPI::broadcast (this->get_mpi_communicator(), true, 0); + // if (openlem_mesh_filled != true) + // throw aspect::QuietException(); + + // // This is called solely so we can set the timer and will return immediately. + // execute_openlem(grid_new, + // //mesh_velocity_z, + // //mesh_velocity_z, + // //mesh_velocity_z, + // //mesh_velocity_z, + // //mesh_velocity_z, + // aspect_timestep_in_years, + // openlem_steps_per_aspect_step); + + //mesh_velocity_z = Utilities::MPI::broadcast(this->get_mpi_communicator(), mesh_velocity_z, 0); + } + //std::cout << "before broadcast mesh_velocity_z" << std::endl; + aspect_mesh_dh = Utilities::MPI::broadcast(this->get_mpi_communicator(), aspect_mesh_dh, 0); + //std::cout << "after broadcast mesh_velocity_z, aspect_dx:dy = " << aspect_dx << ":" << aspect_dy << ", ge = " << grid_extent[0].first << ":" << grid_extent[1].first << ", openlem_nx = " << openlem_nx << ", aspect_nx = " << aspect_nx << std::endl; + //std::cout << "min:max h = " << min_openlem_h << ":" << max_openlem_h << std::endl; + // Get the sizes needed for a data table of the mesh velocities. + //TableIndices size_idx; + //for (unsigned int d=0; d velocity_table;// = fill_data_table(mesh_velocity_z, size_idx, openlem_nx, openlem_ny); + //velocity_table.TableBase::reinit(size_idx); + //// Indexes through x, y, and z. + //if constexpr (dim == 2) + // { + // for (unsigned int x=0; xget_mpi_communicator()) == 0) + connector.write_vtk(grid_extent[dim-1].second, this->get_timestep_number(),this->get_time()/year_in_seconds, dirname); + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) == 1) + connector.write_vtk(grid_extent[dim-1].second, this->get_timestep_number(),this->get_time()/year_in_seconds, dirname, "mpi1"); + //std::cout << velocity_table.size()[0] << ":" << velocity_table.size()[1] << " : " << velocity_table.size()[2] << ", aplha = " << connector.alpha << ", dalpha = " << connector.dalpha << std::endl; + // we get time passed as seconds (always) but may want + // to reinterpret it in years + + // Todo: height_diff = grid_new->height - grid_old->height; + // Todo: interpolate hight_diff to velocties on nodes + // Todo: move info to all other processes + // Todo: copy grid_new to grid_old + } + + + + /** + * A function that creates constraints for the velocity of certain mesh + * vertices (e.g. the surface vertices) for a specific boundary. + * The calling class will respect + * these constraints when computing the new vertex positions. + */ + template + void + OpenLEM::compute_velocity_constraints_on_boundary(const DoFHandler &mesh_deformation_dof_handler, + AffineConstraints &mesh_velocity_constraints, + const std::set &boundary_ids) const + { + if (this->get_timestep_number() == 0) + return; + // Loop over all boundary indicators to set the velocity constraints + /*for (const auto boundary_id : boundary_ids) + VectorTools::interpolate_boundary_values (this->get_mapping(), + mesh_deformation_dof_handler, + boundary_id, + function, + mesh_velocity_constraints); + */ + + + // As our grid_extent variable end points do not account for the change related to an origin + // not at 0, we adjust this here into an interpolation extent. + //std::array,dim> interpolation_extent; + //for (unsigned int d=0; d *velocities; + //Functions::InterpolatedUniformGridData velocities (interpolation_extent, + // table_intervals, + // velocity_table); + + VectorFunctionFromScalarFunctionObject vector_function_object( + [&](const Point &p) -> double + { + //if (p[0] > -6345.87-100. && p[0] < -6345.87+15000. && p[1] > 28103. - 400. && p[1] < 28103. + 400. && p[2] > 25076.5 - 400. && p[2] < 25076.5+5000.) + // //if (p[0] > 3125.00-100. && p[0] < 3125.00+5000. && p[1] > 34375. - 400. && p[1] < 34375. + 400. && p[2] > 28028.4 - 400. && p[2] < 28028.4+5000.) + // { + // //std::cout << "p = " << p << ", vel = " << (connector.interpolation(p[0],p[1], aspect_mesh_elevation_h, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first)-(p[2]-grid_extent[2].second)-(p[2]-grid_extent[2].second))/this->get_timestep() << ", interp = " << connector.interpolation(p[0],p[1], aspect_mesh_elevation_h, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first)-(p[2]-grid_extent[2].second) << ", p2 = " << p[2] << ", ge = " << grid_extent[2].second << ", dt = " << this->get_timestep() << ", without timestep = " << (connector.interpolation(p[0],p[1], aspect_mesh_elevation_h, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first)-(p[2]-grid_extent[2].second)-(p[2]-grid_extent[2].second)) << ", interp/timestep = " << (connector.interpolation(p[0],p[1], aspect_mesh_elevation_h, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first))/this->get_timestep()<< std::endl; + // } + //if ((connector.interpolation(p[0],p[1], aspect_mesh_elevation_h, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first)+ grid_extent[dim-1].second )< 0.0) + //if ((connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first) )> 0.0) + // { + // std::cout << "rank = " << Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) << "p = " << p << ", vel = " << (connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first))/this->get_timestep() << ", interp = " << connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first) << ", p2 = " << p[2] << ", ge = " << grid_extent[2].second << ", dt = " << this->get_timestep() << ", without timestep = " << (connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first)-(p[2]-grid_extent[2].second)-(p[2]-grid_extent[2].second)) << ", interp/timestep = " << (connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first))/this->get_timestep()<< std::endl; + // // AssertThrow(false, ExcMessage("wrong!!!!!!!")) ; + // } + return (connector.interpolation(p[0],p[1], aspect_mesh_dh, aspect_mesh_z, aspect_dx, aspect_dy, grid_extent[0].first, grid_extent[1].first))/this->get_timestep();//+ grid_extent[dim-1].second )/this->get_timestep();//-(p[2]-grid_extent[2].second))/this->get_timestep();// velocities.value(p); + }, + dim-1, + dim); + + VectorTools::interpolate_boundary_values (mesh_deformation_dof_handler, + *boundary_ids.begin(), + vector_function_object, + mesh_velocity_constraints); + + } + + template + std::vector> + OpenLEM::get_aspect_values() const + { + + /* + const types::boundary_id relevant_boundary = this->get_geometry_model().translate_symbolic_boundary_name_to_id ("top"); + std::vector> local_aspect_values(dim+3, std::vector()); + + // Get a quadrature rule that exists only on the corners, and increase the refinement if specified. + const QIterated face_corners (QTrapezoid<1>(), + Utilities::pow(2,additional_refinement_levels+surface_refinement_difference)); + + FEFaceValues fe_face_values (this->get_mapping(), + this->get_fe(), + face_corners, + update_values | + update_quadrature_points); + + for (const auto &cell : this->get_dof_handler().active_cell_iterators()) + if (cell->is_locally_owned() && cell->at_boundary()) + for (unsigned int face_no = 0; face_no < GeometryInfo::faces_per_cell; ++face_no) + if (cell->face(face_no)->at_boundary()) + { + if ( cell->face(face_no)->boundary_id() != relevant_boundary) + continue; + + std::vector> vel(face_corners.size()); + fe_face_values.reinit(cell, face_no); + fe_face_values[this->introspection().extractors.velocities].get_function_values(this->get_solution(), vel); + + for (unsigned int corner = 0; corner < face_corners.size(); ++corner) + { + const Point vertex = fe_face_values.quadrature_point(corner); + + // Find what x point we're at. Add 1 or 2 depending on if ghost nodes are used. + // Subtract the origin point so that it corresponds to an origin of 0,0 in openlem. + const double indx = (vertex(0) - grid_extent[0].first)/openlem_dx; + + // The quadrature rule is created so that there are enough interpolation points in the + // lowest resolved ASPECT surface cell to fill out the openlem mesh. However, as the + // same rule is used for all cell sizes, higher resolution areas will have interpolation + // points that do not correspond to a openlem node. In which case, indx will not be a + // whole number and we can ignore the point. + if (std::abs(indx - std::round(indx)) >= node_tolerance) + continue; + + + // If we're in 2D, we want to take the values and apply them to every row of X points. + if (dim == 2) + { + for (unsigned int ys=0; ys= node_tolerance) + continue; + + //const double index = std::round((indy-1))*openlem_nx+std::round(indx); + + local_aspect_values[0].push_back(vertex(dim-1) - grid_extent[dim-1].second); //z component + if (indx == 5 && indy == 7) + std::cout << "vertex(dim-1) = " << vertex(dim-1) << ", ge[dim-1] = " << grid_extent[dim-1].second << std::endl; + local_aspect_values[1].push_back(indx); + local_aspect_values[2].push_back(indy); + + for (unsigned int d=0; d>(); + } + template + void OpenLEM::execute_openlem(openlem::OceanGrid &grid, + //std::vector &elevation, + //std::vector &extra_vtk_field, + //std::vector &velocity_x, + //std::vector &velocity_y, + //std::vector &velocity_z, + const double &openlem_timestep_in_years, + const unsigned int &openlem_iterations, + const std::string &dirname) + { + TimerOutput::Scope timer_section(this->get_computing_timer(), "Execute openlem"); + if (Utilities::MPI::this_mpi_process(this->get_mpi_communicator()) != 0) + return; + + // Because on the first timestep we will create an initial VTK file before running openlem + // and a second after, we first set the visualization step to zero. + unsigned int visualization_step = 0; + const unsigned int current_timestep = this->get_timestep_number (); + + // Set time step + //openlem_set_dt_(&openlem_timestep_in_years); + this->get_pcout() << " Executing openlem... " << (openlem_iterations) << " timesteps of " << openlem_timestep_in_years << " years." << std::endl; + { + // If it is the first timestep, write an initial VTK file. + if (current_timestep == 1) + { + this->get_pcout() << " Writing initial VTK..." << std::endl; + // openlem by default visualizes a field called HHHHH, + // and the parameter this shows will be whatever is given as the first + // position. At the moment it visualizes the bedrock diffusivity. + //openlem_named_vtk_(extra_vtk_field.data(), + // &vexp, + // &visualization_step, + // dirname_char, + // &dirname_length); + } + constexpr double sea_level = 0; + + for (unsigned int openlem_iteration = 0; openlem_iteration < openlem_iterations; ++openlem_iteration)//openlem_iterations; ++openlem_iteration) + { + //std::cout << "before = " << (grid.getNode(5,7))->h << std::endl; + grid.computeWaterLevel(); + grid.clearOceans(); + grid.markOcean(deepest_point, sea_level); + int nc = grid.computeFlowDirection(); + double ch = grid.erode(openlem_timestep_in_years); + grid.findDeltas(openlem_timestep_in_years); + grid.diffuse(openlem_timestep_in_years); + printf("it: %i, Changes in flow direction: %i, Maximum elevation change: %e; ",openlem_iteration,nc,ch); + unsigned int aspect_timestep_number = this->get_timestep_number(); + connector.updateCoordinateSystem (openlem_timestep_in_years); + //connector.write_vtk(grid_extent[dim-1].second, aspect_timestep_number*10000+openlem_iteration,dirname); + //std::cout << "after = " << (grid.getNode(5,7))->h << std::endl; + //openlem_execute_step_(); + + //// If we are using the ghost nodes we want to reset them every openlem timestep. + //if (use_ghost_nodes) + // { + // openlem_copy_h_(elevation.data()); + + // set_ghost_nodes(elevation, + // velocity_x, + // velocity_y, + // velocity_z, + // openlem_timestep_in_years, + // false); + + // // Set velocity components. + // if (openlem_advection_uplift) + // { + // openlem_set_u_(velocity_z.data()); + // openlem_set_v_(velocity_x.data(), + // velocity_y.data()); + // } + + // // Set h to new values, and erosional parameters if there have been changes. + // openlem_set_h_(elevation.data()); + // } + } + std::cout << "deepest point = " << deepest_point.i << ":" << deepest_point.j << ", deepest point q = " << grid_new[deepest_point].q << ", deepest point h = " << grid_new[deepest_point].h << ", b = " << (int)grid_new[deepest_point].b<< std::endl; + //updateCoordinateSystem + // Copy h values. + //openlem_copy_h_(elevation.data()); + + + // Determine whether to create a VTK file this timestep. + bool write_vtk = false; + + // if (this->get_time() >= last_output_time + output_interval || this->get_time() == this->get_end_time()) + // { + // write_vtk = true; + + // 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; + // } + // } + + // if (write_vtk) + // { + // this->get_pcout() << " Writing openlem VTK..." << std::endl; + // visualization_step = current_timestep; + // openlem_named_vtk_(extra_vtk_field.data(), + // &vexp, + // &visualization_step, + // dirname_char, + // &dirname_length); + // } + } + } + + + template + void OpenLEM::fill_openlem_arrays(openlem::OceanGrid &grid, + // std::vector &elevation, + // std::vector &bedrock_transport_coefficient_array, + // std::vector &bedrock_river_incision_rate_array, + // std::vector &velocity_x, + // std::vector &velocity_y, + // std::vector &velocity_z, + std::vector> &local_aspect_values) + { + for (unsigned int i=0; ih = local_aspect_values[0][i]; + node->l = 0; + node->u = 2000;//1;//local_aspect_values[dim+2][i]; + //if (index_x == 5 && index_y == 7) + // std::cout << "Flag fill 1 before = " << node->h << ", " << local_aspect_values[0][i] << ", node_u = " << node->u << std::endl; + //elevation[index] = local_aspect_values[0][i]; + //velocity_x[index] = local_aspect_values[2][i]; + //velocity_z[index] = local_aspect_values[dim+1][i]; + + //if (dim == 2) + // velocity_y[index] = 0; + //else + // velocity_y[index] = local_aspect_values[3][i]; + } + + for (unsigned int p=1; pget_mpi_communicator()); ++p) + { + // First, find out the size of the array a process wants to send. + MPI_Status status; + MPI_Probe(p, 42, this->get_mpi_communicator(), &status); + int incoming_size = 0; + MPI_Get_count(&status, MPI_DOUBLE, &incoming_size); + + // Resize the array so it fits whatever the process sends. + for (unsigned int i=0; iget_mpi_communicator(), &status); + + // Now, place the numbers into the correct place based off the index. + for (unsigned int i=0; ih = local_aspect_values[0][i]; + node->l = 0; + node->u = 1000;//local_aspect_values[dim+2][i]; + //if (index_x == 5 && index_y == 7) + // std::cout << "Flag fill 2 before = " << node->h << ", " << local_aspect_values[0][i] << ", node_u = " << node->u << std::endl; + //const unsigned int index = local_aspect_values[1][i]; + //elevation[index] = local_aspect_values[0][i]; + //velocity_x[index] = local_aspect_values[2][i]; + //velocity_z[index] = local_aspect_values[dim+1][i]; + + //// In 2D there are no y velocities, so we set them to zero. + //if (dim == 2 ) + // velocity_y[index] = 0; + //else + // velocity_y[index] = local_aspect_values[3][i]; + } + } + + // Initialize the bedrock river incision rate and transport coefficient, + // and check that there are no empty mesh points due to + // an improperly set maximum_surface_refinement_level, additional_refinement_levels, + // and surface_refinement_difference + bool openlem_mesh_filled = true; + const unsigned int openlem_array_size = openlem_nx*openlem_ny; + for (unsigned int i=0; i::max() && !is_ghost_node(i,false)) + // openlem_mesh_filled = false; + } + + Utilities::MPI::broadcast(this->get_mpi_communicator(), openlem_mesh_filled, 0); + AssertThrow (openlem_mesh_filled == true, + ExcMessage("The openlem mesh is missing data. A likely cause for this is that the " + "maximum surface refinement or surface refinement difference are improperly set.")); + } + + + + template + bool + OpenLEM:: + needs_surface_stabilization () const + { + return false; + } + + + + template + void OpenLEM::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection ("Mesh deformation"); + { + prm.enter_subsection ("OpenLEM"); + { + prm.declare_entry("Minimize advection", "true", + Patterns::Bool(), + "Whether to rotate and translate the openLEM grid to minimize the advection"); + prm.declare_entry("Number of openlem timesteps per aspect timestep", "5", + Patterns::Integer(), + "Initial number of openlem time steps per ASPECT timestep, this value will double if" + "the openlem timestep is above the maximum openlem timestep."); + prm.declare_entry("Kd", "1e-14", + Patterns::Double(0), + "Kd erosion rate for openlem. Units: $\\{yrs}$"); + prm.declare_entry("Kt", "1e-14", + Patterns::Double(0), + "Kt erosion rate for openlem. Units: $\\{yrs}$"); + prm.declare_entry("Ocean diffusivity", "1", + Patterns::Double(0), + "The diffusivisty of the ocean openlem. Units: $\\{yrs}$"); + prm.declare_entry("Maximum timestep length", "10e3", + Patterns::Double(0), + "Maximum timestep for openlem. Units: $\\{yrs}$"); + prm.declare_entry("Vertical exaggeration", "-1", + Patterns::Double(), + "Vertical exaggeration for openlem's VTK file. -1 outputs topography, basement, and sealevel."); + prm.declare_entry("Additional openlem refinement", "0", + Patterns::Integer(), + "How many levels above ASPECT openlem should be refined."); + prm.declare_entry ("Average out of plane surface topography in 2d", "true", + Patterns::Bool (), + "If this is set to false, then a 2D model will only consider the " + "center slice openlem gives. If set to true, then ASPECT will" + "average the mesh along Y excluding the ghost nodes."); + prm.declare_entry("openlem seed", "1000", + Patterns::Integer(), + "Seed used for adding an initial noise to openlem topography based on the initial noise magnitude."); + prm.declare_entry("Maximum surface refinement level", "1", + Patterns::Integer(), + "This should be set to the highest ASPECT refinement level expected at the surface."); + prm.declare_entry("Surface refinement difference", "0", + Patterns::Integer(), + "The difference between the lowest and highest refinement level at the surface. E.g., if three resolution " + "levels are expected, this would be set to 2."); + prm.declare_entry ("Use marine component", "false", + Patterns::Bool (), + "Flag to use the marine component of openlem."); + prm.declare_entry("Y extent in 2d", "100000", + Patterns::Double(), + "openlem Y extent when using a 2D ASPECT model. Units: $\\{m}$"); + prm.declare_entry ("Use ghost nodes", "true", + Patterns::Bool (), + "Flag to use ghost nodes"); + prm.declare_entry ("Uplift and advect with openlem", "true", + Patterns::Bool (), + "Flag to use openlem advection and uplift."); + prm.declare_entry("Node tolerance", "0.001", + Patterns::Double(), + "Node tolerance for how close an ASPECT node must be to the openlem node for the value to be transferred."); + prm.declare_entry ("Sediment rain rates", "0,0", + Patterns::List (Patterns::Double(0)), + "Sediment rain rates given as a list 1 greater than the number of sediment rain time intervals. E.g, " + " If the time interval is given at 5 Myr, there will be one value for 0-5 Myr model time and a second value " + " for 5+ Myr. Units: $\\{m/yr}$"); + prm.declare_entry ("Sediment rain time intervals", "0", + Patterns::List (Patterns::Double(0)), + "A list of times to change the sediment rain rate. Units: $\\{yrs}$"); + prm.declare_entry("Initial noise magnitude", "5", + Patterns::Double(), + "Maximum topography change from the initial noise. Units: $\\{m}$"); + + prm.enter_subsection ("Boundary conditions"); + { + prm.declare_entry ("Front", "1", + Patterns::Integer (0, 1), + "Front (bottom) boundary condition, where 1 is fixed and 0 is reflective."); + prm.declare_entry ("R ight", "1", + Patterns::Integer (0, 1), + "Right boundary condition, where 1 is fixed and 0 is reflective."); + prm.declare_entry ("Back", "1", + Patterns::Integer (0, 1), + "Back (top) boundary condition, where 1 is fixed and 0 is reflective."); + prm.declare_entry ("Left", "1", + Patterns::Integer (0, 1), + "Left boundary condition, where 1 is fixed and 0 is reflective."); + prm.declare_entry("Left mass flux", "0", + Patterns::Double(), + "Flux per unit length through left boundary. Units: $\\{m^2/yr}$ "); + prm.declare_entry("Right mass flux", "0", + Patterns::Double(), + "Flux per unit length through right boundary. Units: $\\{m^2/yr}$ "); + prm.declare_entry("Back mass flux", "0", + Patterns::Double(), + "Flux per unit length through back boundary. Units: $\\{m^2/yr}$ "); + prm.declare_entry("Front mass flux", "0", + Patterns::Double(), + "Flux per unit length through front boundary. Units: $\\{m^2/yr}$ "); + prm.declare_entry ("Back front ghost nodes periodic", "false", + Patterns::Bool (), + "Whether to set the ghost nodes at the openlem back and front boundary " + "to periodic even if 'Back' and 'Front' are set to fixed boundary."); + prm.declare_entry ("Left right ghost nodes periodic", "false", + Patterns::Bool (), + "Whether to set the ghost nodes at the openlem left and right boundary " + "to periodic even if 'Left' and 'Right' are set to fixed boundary."); + } + prm.leave_subsection(); + + prm.enter_subsection ("Erosional parameters"); + { + prm.declare_entry("Drainage area exponent", "0.4", + Patterns::Double(), + "Exponent for drainage area."); + prm.declare_entry("Slope exponent", "1", + Patterns::Double(), + "The slope exponent for SPL (n). Generally m/n should equal approximately 0.4"); + prm.declare_entry("Multi-direction slope exponent", "1", + Patterns::Double(), + "Exponent to determine the distribution from the SPL to neighbor nodes, with" + "10 being steepest decent and 1 being more varied."); + prm.declare_entry("Bedrock deposition coefficient", "1", + Patterns::Double(), + "Deposition coefficient for bedrock."); + prm.declare_entry("Sediment deposition coefficient", "-1", + Patterns::Double(), + "Deposition coefficient for sediment, -1 sets this to the same as the bedrock deposition coefficient."); + prm.declare_entry("Bedrock river incision rate", "1e-5", + Patterns::Double(), + "River incision rate for bedrock in the Stream Power Law. Units: $\\{m^(1-2*drainage_area_exponent)/yr}$"); + prm.declare_entry("Sediment river incision rate", "-1", + Patterns::Double(), + "River incision rate for sediment in the Stream Power Law. -1 sets this to the bedrock river incision rate. Units: $\\{m^(1-2*drainage_area_exponent)/yr}$ "); + prm.declare_entry("Bedrock diffusivity", "1e-2", + Patterns::Double(), + "Transport coefficient (diffusivity) for bedrock. Units: $\\{m^2/yr}$ "); + prm.declare_entry("Sediment diffusivity", "-1", + Patterns::Double(), + "Transport coefficient (diffusivity) for sediment. -1 sets this to the bedrock diffusivity. Units: $\\{m^2/yr}$"); + prm.declare_entry("Orographic elevation control", "2000", + Patterns::Integer(), + "Above this height, the elevation factor is applied. Units: $\\{m}$"); + prm.declare_entry("Orographic wind barrier height", "500", + Patterns::Integer(), + "When terrain reaches this height the wind barrier factor is applied. Units: $\\{m}$"); + prm.declare_entry("Elevation factor", "1", + Patterns::Double(), + "Amount to multiply the bedrock river incision rate nad transport coefficient by past the given orographic elevation control."); + prm.declare_entry("Wind barrier factor", "1", + Patterns::Double(), + "Amount to multiply the bedrock river incision rate nad transport coefficient by past given wind barrier height."); + prm.declare_entry ("Stack orographic controls", "true", + 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 ("Flag to use orographic controls", "false", + Patterns::Bool (), + "Whether or not to apply orographic controls."); + 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."); + prm.declare_entry ("Use a fixed erosional base 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 openlem boundaries to the height you specify in " + "'set Erosional base 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 openlem 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 openlem."); + prm.declare_entry("Erosional base level", "0", + Patterns::Double(), + "When 'Use a fixed erosional base level' is set to true, all ghost nodes of fixed " + "openlem boundaries where no mass flux is specified by the user (openlem boundary condition set to 1 " + "and 'Left/Right/Bottom/Top mass flux' set to 0) will be fixed to this elevation. The " + "reflecting boundaries (openlem 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(); + + prm.enter_subsection ("Marine parameters"); + { + prm.declare_entry("Sea level", "0", + Patterns::Double(), + "Sea level relative to the ASPECT surface, where the maximum Z or Y extent in ASPECT is a sea level of zero. Units: $\\{m}$ "); + prm.declare_entry("Sand porosity", "0.0", + Patterns::Double(), + "Porosity of sand. "); + prm.declare_entry("Silt porosity", "0.0", + Patterns::Double(), + "Porosity of silt. "); + prm.declare_entry("Sand e-folding depth", "1e3", + Patterns::Double(), + "E-folding depth for the exponential of the sand porosity law. Units: $\\{m}$"); + prm.declare_entry("Silt e-folding depth", "1e3", + Patterns::Double(), + "E-folding depth for the exponential of the silt porosity law. Units: $\\{m}$"); + prm.declare_entry("Sand-silt ratio", "0.5", + Patterns::Double(), + "Ratio of sand to silt for material leaving continent."); + prm.declare_entry("Depth averaging thickness", "1e2", + Patterns::Double(), + "Depth averaging for the sand-silt equation. Units: $\\{m}$"); + prm.declare_entry("Sand transport coefficient", "5e2", + Patterns::Double(), + "Transport coefficient (diffusivity) for sand. Units: $\\{m^2/yr}$"); + prm.declare_entry("Silt transport coefficient", "2.5e2", + Patterns::Double(), + "Transport coefficient (diffusivity) for silt. Units: $\\{m^2/yr}$ "); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection (); + } + + template + void OpenLEM::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection ("Mesh deformation"); + { + prm.enter_subsection("OpenLEM"); + { + openlem_ocean_diffusivity = prm.get_double("Ocean diffusivity"); + openlem_minimize_advection = prm.get_bool("Minimize advection"); + openlem_kd = prm.get_double("Kd"); + openlem_kt = prm.get_double("Kt"); + openlem_steps_per_aspect_step = prm.get_integer("Number of openlem timesteps per aspect timestep"); + maximum_openlem_timestep = prm.get_double("Maximum timestep length"); + //vexp = prm.get_double("Vertical exaggeration"); + additional_refinement_levels = prm.get_integer("Additional openlem refinement"); + //average_out_of_plane_surface_topography = prm.get_bool("Average out of plane surface topography in 2d"); + //openlem_seed = prm.get_integer("openlem seed"); + maximum_surface_refinement_level = prm.get_integer("Maximum surface refinement level"); + surface_refinement_difference = prm.get_integer("Surface refinement difference"); + //use_marine_component = prm.get_bool("Use marine component"); + openlem_y_extent_2d = prm.get_double("Y extent in 2d"); + use_ghost_nodes = prm.get_bool("Use ghost nodes"); + //openlem_advection_uplift = prm.get_bool("Uplift and advect with openlem"); + node_tolerance = prm.get_double("Node tolerance"); + //noise_elevation = prm.get_double("Initial noise magnitude"); + //sediment_rain_rates = Utilities::string_to_double + // (Utilities::split_string_list(prm.get ("Sediment rain rates"))); + //sediment_rain_times = Utilities::string_to_double + // (Utilities::split_string_list(prm.get ("Sediment rain time intervals"))); + + //if (!this->convert_output_to_years()) + // { + // maximum_openlem_timestep /= year_in_seconds; + // for (unsigned int j=0; j sediment_rain_times[i-1], ExcMessage("Sediment rain time intervals must be an increasing array.")); + + prm.enter_subsection("Boundary conditions"); + { + // bottom = prm.get_integer("Front"); + // right = prm.get_integer("Right"); + // top = prm.get_integer("Back"); + // left = prm.get_integer("Left"); + // left_flux = prm.get_double("Left mass flux"); + // right_flux = prm.get_double("Right mass flux"); + // top_flux = prm.get_double("Back mass flux"); + // bottom_flux = prm.get_double("Front mass flux"); + + // if (!this->convert_output_to_years()) + // { + // left_flux *= year_in_seconds; + // right_flux *= year_in_seconds; + // top_flux *= year_in_seconds; + // bottom_flux *= year_in_seconds; + // } + + // // Put the boundary condition values into a four digit value to send to openlem. + // openlem_boundary_conditions = bottom*1000+right*100+top*10+left; + + // if ((left_flux != 0 && top_flux != 0) || (left_flux != 0 && bottom_flux != 0) || + // (right_flux != 0 && bottom_flux != 0) || (right_flux != 0 && top_flux != 0)) + // AssertThrow(false,ExcMessage("Currently the plugin does not support mass flux through adjacent boundaries.")); + + // topbottom_ghost_nodes_periodic = prm.get_bool("Back front ghost nodes periodic"); + // leftright_ghost_nodes_periodic = prm.get_bool("Left right ghost nodes periodic"); + } + prm.leave_subsection(); + + prm.enter_subsection("Erosional parameters"); + { + //drainage_area_exponent_m = prm.get_double("Drainage area exponent"); + //slope_exponent_n = prm.get_double("Slope exponent"); + //sediment_river_incision_rate = prm.get_double("Sediment river incision rate"); + //bedrock_river_incision_rate = prm.get_double("Bedrock river incision rate"); + //sediment_transport_coefficient = prm.get_double("Sediment diffusivity"); + //bedrock_transport_coefficient = prm.get_double("Bedrock diffusivity"); + //bedrock_deposition_g = prm.get_double("Bedrock deposition coefficient"); + //sediment_deposition_g = prm.get_double("Sediment deposition coefficient"); + //slope_exponent_p = prm.get_double("Multi-direction slope exponent"); + //flat_elevation = prm.get_integer("Orographic elevation control"); + //wind_barrier_elevation = prm.get_integer("Orographic wind barrier height"); + //flat_erosional_factor = prm.get_double("Elevation factor"); + //wind_barrier_erosional_factor = prm.get_double("Wind barrier factor"); + //stack_controls = prm.get_bool("Stack orographic controls"); + //use_orographic_controls = prm.get_bool("Flag to use orographic controls"); + + //if (!this->convert_output_to_years()) + // { + // bedrock_river_incision_rate *= year_in_seconds; + // bedrock_transport_coefficient *= year_in_seconds; + // sediment_river_incision_rate *= year_in_seconds; + // bedrock_transport_coefficient *= year_in_seconds; + // } + + //// Wind direction + //if (prm.get ("Wind direction") == "west") + // wind_direction = 0; + //else if (prm.get ("Wind direction") == "east") + // wind_direction = 1; + //else if (prm.get ("Wind direction") == "north") + // wind_direction = 2; + //else if (prm.get ("Wind direction") == "south") + // wind_direction = 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_fixed_erosional_base = prm.get_bool("Use a fixed erosional base level"); + //if (use_fixed_erosional_base) + // AssertThrow(use_fixed_erosional_base && use_ghost_nodes, ExcMessage( + // "If you want to use an erosional base level differing from sea level, " + // "you need to use ghost nodes.")); + //h_erosional_base = prm.get_double("Erosional base level"); + } + prm.leave_subsection(); + + prm.enter_subsection("Marine parameters"); + { + //sea_level = prm.get_double("Sea level"); + //sand_surface_porosity = prm.get_double("Sand porosity"); + //silt_surface_porosity = prm.get_double("Silt porosity"); + //sand_efold_depth = prm.get_double("Sand e-folding depth"); + //silt_efold_depth = prm.get_double("Silt e-folding depth"); + //sand_silt_ratio = prm.get_double("Sand-silt ratio"); + //sand_silt_averaging_depth = prm.get_double("Depth averaging thickness"); + //sand_transport_coefficient = prm.get_double("Sand transport coefficient"); + //silt_transport_coefficient = prm.get_double("Silt transport coefficient"); + + //if (!this->convert_output_to_years()) + // { + // sand_transport_coefficient *= year_in_seconds; + // silt_transport_coefficient *= year_in_seconds; + // } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + prm.leave_subsection (); + } + } +} + + +// explicit instantiation of the functions we implement in this file +namespace aspect +{ + namespace MeshDeformation + { + ASPECT_REGISTER_MESH_DEFORMATION_MODEL(OpenLEM, + "openLEM", + "TODO: A plugin, which prescribes the surface mesh to " + "deform according to an analytically prescribed " + "function. Note that the function prescribes a " + "deformation velocity, i.e. the return value of " + "this plugin is later multiplied by the time step length " + "to compute the displacement increment in this time step. " + "Although the function's time variable is interpreted as " + "years when Use years in output instead of seconds is set to true, " + "the boundary deformation velocity should still be given " + "in m/s. The format of the " + "functions follows the syntax understood by the " + "muparser library, see {ref}\\`sec:run-aspect:parameters-overview:muparser-format\\`.") + } +}