Prozedurale Hydrologie: Dynamische Simulation von Flüssen und Seen

Hinweis: Der vollständige Quellcode des Projekts wird auf Github [ hier ] veröffentlicht. Das Repository enthält außerdem detaillierte Informationen zum Lesen und Verwenden des Codes.

Erläuterung: Ich bin kein Geologe, daher habe ich das System basierend auf meinem Wissen erstellt.

Hydrologiekonzept

• Fluss- und Bachwanderung
• Natürliche Wasserfälle
• Canyon Formation
• Schwellung des Bodens und der Auen

• Das Relief beeinflusst die Bewegung des Wassers
• Erosion und Sedimentation beeinflussen das Gelände

• In einem sich bewegenden Strom verhält sich Wasser anders.
• In einem stehenden Pool verhält sich das Wasser anders

Vereinfachtes hydrologisches Modell

Hinweis: Vergessen Sie nicht, dass sie als 2D-Systeme modelliert sind.

Stream- und Pool-Karten

Hinweis: Stream- und Pool-Maps sind Arrays mit derselben Größe wie die Höhenkarte.

Wasser als Partikel

Hydrologische Zyklus- und Karteninteraktion

• Erstellen eines Partikels im Gelände
• (.. ).
• , .
• , .
• , ( ) .
• .

Implementierung

Partikelklasse

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);
};

Abstiegsalgorithmus

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]);
//...

Hinweis: Ich habe herausgefunden, dass eine solche Parameteränderung gut funktioniert.

//... 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;

//...

Hinweis: Durch die Überlaufbedingung wird auch das Volumen zurückgesetzt, sodass die Partikel nicht zum Flutungsalgorithmus übergehen.

Hochwasseralgorithmus

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);

    //...

Hinweis: Der Einfachheit halber wird hier ein Acht- Wege-Füllalgorithmus verwendet . In Zukunft wird es möglich sein, es effizienter umzusetzen.

    //...

    //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;
    }

    //...

Hinweis: Wenn der Wasserstand aufgrund von Leckagen abnimmt, habe ich festgestellt, dass dies bei einer geringen Leckrate am besten funktioniert. Darüber hinaus hilft die Entfernung eines Teils des Sedimentgesteins bei der Umsetzung.

    //...

    //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

Erosionsfolie

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};

};

Hinweis: Die Höhenkarte wird mit Perlin-Rauschen initialisiert.

//...

//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--;
}

//...

Hinweis: Partikel gelangen selten in die Pools und verlassen sie mehr als ein oder zwei Mal, bevor sie während der Abstiegsphasen vollständig verdunsten. Ich habe dies jedoch nur für den Fall hinzugefügt.

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);

}

Bäume

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

Hinweis: Das Baummodell befindet sich in der Datei "vegetation.h" und in der Funktion "World :: grow ()".

Andere Details

Ergebnisse

Vergleich der Auswirkungen der Verengung von Strömungen

• Partikelbasierte Erosion (Grunderosion), die Fluss- und Poolkarten erhält. Pools wirken sich immer noch auf die Generierung aus
• Grundlegende Erosion mit Parametern, die durch hydrologische Karten geändert wurden (in Kombination mit Erosion)
• Kombinierte Erosion mit Parametern, die durch hydrologische Karten geändert wurden und die Erosion durch Bäume beeinflussen

Auswirkungen auf den Pool

Hinweis: Ich habe mehrere Samen erhalten, in denen drei Seen nacheinander ineinander flossen, aber ich wollte nicht viel Zeit damit verbringen, den richtigen für diesen Artikel zu finden. Ich habe bereits zu viele Bilder und Videos generiert.

Simulationsgeschwindigkeit

Hinweis: Kurz nach der Veröffentlichung dieses Artikels habe ich die Methode zur Berechnung von Oberflächennormalen geringfügig geändert, wodurch die Simulationsgeschwindigkeit erhöht wurde. Der Effekt ist während der Simulation sehr deutlich, aber es ist schwierig, ihn zu bewerten, da die Zeit zum Starten und Fluten sehr unterschiedlich ist. Nach meinen Schätzungen hat sich die Simulationsgeschwindigkeit verdoppelt.

Weitere schöne Videos

Fazit und Arbeit für die Zukunft