Hidrologia processual: simulação dinâmica de rios e lagos

Nota: o código fonte completo do projeto está publicado no Github [ aqui ]. O repositório também contém informações detalhadas sobre como ler e usar o código.

Explicação: Como não sou geólogo, criei o sistema com base no meu conhecimento.

Conceito de hidrologia

• Migração de rios e riachos
• Cachoeiras naturais
• Formação Canyon
• Inchaço do solo e várzea

• O alívio afeta o movimento da água
• Erosão e sedimentação afetam o terreno

• Em um fluxo em movimento , a água se comporta de maneira diferente.
• Em uma piscina em pé , a água se comporta de maneira diferente

Modelo hidrológico simplificado

Nota: não esqueça que eles são modelados como sistemas 2D.

Mapas de Stream e Pool

Nota: os mapas de fluxo e pool são matrizes do mesmo tamanho que o mapa de altura.

Água como uma partícula

Ciclo hidrológico e interação do mapa

• Criando uma partícula no terreno
• (.. ).
• , .
• , .
• , ( ) .
• .

Implementação

Classe de partículas

struct Drop{
  
  //... constructors

  int index;                         //Flat Array Index
  glm::vec2 pos;                     //2D Position
  glm::vec2 speed = glm::vec2(0.0);
  double volume = 1.0;
  double sediment = 0.0;

  //... parameters
  const double volumeFactor = 100.0; //"Water Deposition Rate"

  //Hydrological Cycle Functions
  void descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale);
  void flood(double* h, double* pool, glm::ivec2 dim);
};

Algoritmo de descida

void Drop::descend(double* h, double* stream, double* pool, bool* track, glm::ivec2 dim, double scale){

  glm::ivec2 ipos; 

  while(volume > minVol){

    ipos = pos; //Initial Position
    int ind = ipos.x*dim.y+ipos.y; //Flat Array Index

    //Register Position
    track[ind] = true;

    //...
  }
};

//...
  //Effective Parameter Set
  double effD = depositionRate;
  double effF = friction*(1.0-0.5*stream[ind]);
  double effR = evapRate*(1.0-0.2*stream[ind]);
//...

Nota: Descobri que essa alteração de parâmetro funciona bem.

//... nind is the next position after moving the particle
  
  //Out-Of-Bounds
  if(!glm::all(glm::greaterThanEqual(pos, glm::vec2(0))) ||
     !glm::all(glm::lessThan((glm::ivec2)pos, dim))){
       volume = 0.0;
       break;
  }

  //Slow-Down
  if(stream[nind] > 0.5 && length(acc) < 0.01)
    break;

  //Enter Pool
  if(pool[nind] > 0.0)
    break;

//...

Nota: a condição de transbordamento também redefine o volume para que as partículas não passem para o algoritmo de inundação.

Algoritmo de inundação

void Drop::flood(double* height, double* pool, glm::ivec2 dim){

  index = (int)pos.x*dim.y + (int)pos.y;
  double plane = height[index] + pool[index];  //Testing Plane
  double initialplane = plane;                 //Water Level

  //Flood Set
  std::vector<int> set;
  int fail = 10; //Just in case...

  //Iterate while particle still has volume
  while(volume > minVol && fail){

    set.clear();
    bool tried[dim.x*dim.y] = {false};

    //Lowest Drain
    int drain;
    bool drainfound = false;

    //Recursive Flood-Fill Function
    std::function<void(int)> fill = [&](int i){

      //Out of Bounds
      if(i/dim.y >= dim.x || i/dim.y < 0) return;
      if(i%dim.y >= dim.y || i%dim.y < 0) return;

      //Position has been tried
      if(tried[i]) return;
      tried[i] = true;

      //Wall / Boundary of the Pool
      if(plane < height[i] + pool[i]) return;

      //Drainage Point
      if(initialplane > height[i] + pool[i]){

        //No Drain yet
        if(!drainfound)
          drain = i;

        //Lower Drain
        else if( pool[drain] + height[drain] < pool[i] + height[i] )
          drain = i;

        drainfound = true;
        return; //No need to flood from here
      }

      //Part of the Pool
      set.push_back(i);
      fill(i+dim.y);    //Fill Neighbors
      fill(i-dim.y);
      fill(i+1);
      fill(i-1);
      fill(i+dim.y+1);  //Diagonals (Improves Drainage)
      fill(i-dim.y-1);
      fill(i+dim.y-1);
      fill(i-dim.y+1);
    };

    //Perform Flood
    fill(index);

    //...

Nota: para simplificar, um algoritmo de preenchimento de oito vias é usado aqui . No futuro, será possível implementá-lo com mais eficiência.

    //...

    //Drainage Point
    if(drainfound){

      //Set the Particle Position
      pos = glm::vec2(drain/dim.y, drain%dim.y);

      //Set the New Waterlevel (Slowly)
      double drainage = 0.001;
      plane = (1.0-drainage)*initialplane + drainage*(height[drain] + pool[drain]);

      //Compute the New Height
      for(auto& s: set) //Iterate over Set
        pool[s] = (plane > height[s])?(plane-height[s]):0.0;

      //Remove some sediment
      sediment *= 0.1;
      break;
    }

    //...

Nota: quando o nível da água diminui devido a vazamentos, descobri que isso funciona melhor com uma baixa taxa de vazamentos. Além disso, a eliminação de parte das rochas sedimentares ajuda na implementação.

    //...

    //Get Volume under Plane
    double tVol = 0.0;
    for(auto& s: set)
      tVol += volumeFactor*(plane - (height[s]+pool[s]));

    //We can partially fill this volume
    if(tVol <= volume && initialplane < plane){

      //Raise water level to plane height
      for(auto& s: set)
        pool[s] = plane - height[s];

      //Adjust Drop Volume
      volume -= tVol;
      tVol = 0.0;
    }

    //Plane was too high and we couldn't fill it
    else fail--;

    //Adjust Planes
    float approach = 0.5;
    initialplane = (plane > initialplane)?plane:initialplane;
    plane += approach*(volume-tVol)/(double)set.size()/volumeFactor;
  }

  //Couldn't place the volume (for some reason)- so ignore this drop.
  if(fail == 0)
    volume = 0.0;

} //End of Flood Algorithm

Envoltório de erosão

class World {

public:
  void generate();            //Initialize
  void erode(int cycles);     //Erode with N Particles

  //...

  double heightmap[256*256] = {0.0};
  double waterstream[256*256] = {0.0};
  double waterpool[256*256] = {0.0};

};

Nota: O mapa de altura é inicializado usando o ruído Perlin.

//...

//Spawn Particle
glm::vec2 newpos = glm::vec2(rand()%(int)dim.x, rand()%(int)dim.y);
Drop drop(newpos);

int spill = 5;
while(drop.volume > drop.minVol && spill != 0){

  drop.descend(heightmap, waterstream, waterpool, track, dim, scale);

  if(drop.volume > drop.minVol)
    drop.flood(heightmap, waterpool, dim);

  spill--;
}

//...

Nota: as partículas raramente entram nas piscinas e as deixam mais de uma ou duas vezes antes de evaporar completamente durante as etapas de descida, mas adicionei isso apenas por precaução.

void World::erode(int N){

  //Track the Movement of all Particles
  bool track[dim.x*dim.y] = {false};

  //Simulate N Particles
  for(int i = 0; i < N; i++){
   
    //... simulate individual particle

  }

  //Update Path
  double lrate = 0.01;  //Adaptation Rate
  for(int i = 0; i < dim.x*dim.y; i++)
    waterstream[i] = (1.0-lrate)*waterstream[i] + lrate*((track[i])?1.0:0.0);

}

Árvores

//... descend function
double effD = depositionRate*max(0.0, 1.0-treedensity[ind]);
//...

Nota: o modelo da árvore pode ser encontrado no arquivo "vegetação.h" e na função "Mundo :: crescer ()".

Outros detalhes

resultados

Comparação dos efeitos do estreitamento dos fluxos

• Erosão baseada em partículas (erosão básica) que recebe mapas de fluxo e piscina. Piscinas ainda afetam geração
• Erosão básica com parâmetros alterados por mapas hidrológicos (em combinação com erosão)
• Erosão combinada com parâmetros alterados por mapas hidrológicos e afetando a erosão por árvores

Impacto na piscina

Nota: Eu recebi várias sementes, nas quais três lagos corriam um para o outro sequencialmente, mas não queria gastar muito tempo encontrando o caminho certo para este artigo. Eu já gerei muitas fotos e vídeos.

Velocidade de simulação

Nota: logo após a publicação deste artigo, fiz uma pequena alteração no método de cálculo das normais da superfície, o que aumentou a velocidade da simulação. O efeito disso é muito perceptível durante a simulação, mas é difícil compará-lo devido a grandes variações no tempo necessário para iniciar e inundar. De acordo com minhas estimativas, a velocidade da simulação dobrou.

Vídeos mais bonitos

Conclusão e trabalho para o futuro