الهيدرولوجيا الإجرائية: محاكاة ديناميكية للأنهار والبحيرات

ملاحظة: تم نشر كود المصدر الكامل للمشروع على Github [ هنا ]. يحتوي المستودع أيضًا على معلومات تفصيلية حول كيفية قراءة التعليمات البرمجية واستخدامها.

التفسير: لست جيولوجيًا ، لذلك أنشأت النظام بناءً على معرفتي.

مفهوم الهيدرولوجيا

• هجرة الأنهار والجداول
• الشلالات الطبيعية
• تشكيل الوادي
• تورم التربة والسهول الفيضية

• يؤثر التخفيف على حركة الماء
• يؤثر التعرية والترسيب على التضاريس

• في تيار متحرك ، يتصرف الماء بشكل مختلف.
• يتصرف الماء بشكل مختلف في حمام السباحة القائم

نموذج هيدرولوجي مبسط

ملحوظة: لا تنس أنهم مصممون على هيئة أنظمة ثنائية الأبعاد.

دفق وخرائط البلياردو

ملحوظة: خرائط التدفق والتجمع هي صفائف بنفس حجم خريطة الارتفاع.

الماء كجسيم

الدورة الهيدرولوجية وتفاعل الخريطة

• خلق جسيم على التضاريس
• (.. ).
• , .
• , .
• , ( ) .
• .

التنفيذ

فئة الجسيمات

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

خوارزمية النسب

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

ملاحظة: اكتشفت أن تغيير المعلمة يعمل بشكل جيد.

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

//...

ملاحظة: تقوم حالة الفائض أيضًا بإعادة تعيين الحجم بحيث لا تنتقل الجزيئات إلى خوارزمية الفيضان.

خوارزمية الفيضانات

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

    //...

ملاحظة: من أجل البساطة ، يتم استخدام خوارزمية تعبئة ذات ثمانية اتجاهات هنا . في المستقبل ، سيكون من الممكن تنفيذها بشكل أكثر كفاءة.

    //...

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

    //...

ملاحظة: عندما ينخفض ​​مستوى المياه بسبب التسرب ، اكتشفت أن هذا يعمل بشكل أفضل مع معدل تسرب منخفض. بالإضافة إلى ذلك ، يساعد التخلص من جزء من الصخور الرسوبية على التنفيذ.

    //...

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

التفاف التآكل

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

};

ملاحظة: تتم تهيئة خريطة الارتفاع باستخدام ضوضاء 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--;
}

//...

ملاحظة: نادرًا ما تدخل الجسيمات إلى البرك وتتركها أكثر من مرة أو مرتين قبل أن تتبخر تمامًا خلال مراحل الهبوط ، لكني أضفت هذا فقط في حالة.

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

}

الأشجار

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

ملاحظة: يمكن العثور على نموذج الشجرة في الملف "vegetation.h" وفي الوظيفة "World :: Grow ()".

تفاصيل أخرى

النتائج

مقارنة آثار تضييق التدفقات

• التآكل القائم على الجسيمات (التآكل الأساسي) الذي يتلقى خرائط التدفق والتجمع. لا تزال التجمعات تؤثر على التوليد
• تغيرت التآكل الأساسي مع المعلمات بواسطة الخرائط الهيدرولوجية (بالاشتراك مع التآكل)
• التآكل المشترك مع البارامترات التي تغيرتها الخرائط الهيدرولوجية وتؤثر على التآكل بواسطة الأشجار

تأثير التجمع

ملحوظة: حصلت على عدة بذور ، تدفقت فيها ثلاث بحيرات على بعضها بالتسلسل ، لكنني لم أرغب في قضاء الكثير من الوقت في العثور على البذور المناسبة لهذه المقالة. لقد ولدت بالفعل العديد من الصور ومقاطع الفيديو.

سرعة المحاكاة

ملاحظة: بعد فترة وجيزة من نشر هذه المقالة ، قمت بتغيير طفيف في طريقة حساب المعايير السطحية ، مما زاد من سرعة المحاكاة. تأثير هذا ملحوظ للغاية أثناء المحاكاة ، ولكن من الصعب قياسه بسبب الاختلافات الكبيرة في الوقت الذي يستغرقه الإطلاق والفيضان. وفقا لتقديراتي ، تضاعفت سرعة المحاكاة.

فيديوهات أجمل

الخلاصة والعمل من أجل المستقبل