Python rechnen mit Klassen? Objektorientierung

[Tenferenzu]

Liebes Forum,
ich stehe mal wieder auf dem Schlauch was Programmierung angeht. Rein Ablauforientiert bekomm ich das mittlerweile schon ganz gut hin, allerdings blicke ich noch nicht ganz durch wie man wirklich mit Klassen/Instanzen (das sollte der richtige Terminus sein oder?) rechnet.

Es geht hier um einen numerischen Löser der erstmal ganz simpel gehalten werden soll für das Beispiel. Im Prinzip möchte ich jeder Zelle auf meinem (1D) Rechengrid zu einer Instanz machen und ihnen neben einem Index auch ein paar Werte wie Geschwindigkeit u, Temperatur t etc. zuweisen und dann damit rechnen können. Alle Beispiele ich ich so zur Objektorientierung im Bereich Python gefunden habe, waren meist auf andere Anwendungen zugeschnitten. Deshalb hoffe ich hier etwas Hilfe zu bekommen.

Ich nehme einfach mal das Beispiel Geschwindigkeit da das am einfachsten zu berechnen ist. Aktuell habe ich eine Liste/Vektor mit Werten die dann durchiteriert wird um den nachfolgenden Zeitpunkt t=/=0 zu berechnen. Dazu verwende ich die eigentliche Zelle (i) und die danach bzw rechts davon (i+1) als auch zwei Listen. Eine aktuelle Liste und eine für den neuen Zeitpunkt t+1.

Hier mal der Code den ich Beispielhaft überführen möchte:

def linearconv(nt, dt, c, nu, filename, varname):
    u = csvimporttonumpy('%s' %varname)
    un = numpy.ones(nx)
    iterations = [] # neu laden der Variable um alte Werte loszuwerden
    cour = []
    diffzahl = []
    for i in dx:
        courrant = c * dt / i
        cour.append(courrant)
        difz = nu * dt / (i ** 2)
        diffzahl.append(difz)
    print('Min. Courant number is:', min(cour))
    print('Max. Courant number is:', max(cour))
    print('Min. diffusion number is:', min(diffzahl))
    print('Max. diffusion number is:', max(diffzahl))
    for n in range(nt):  # loop for values of n from 0 to nt, so it will run nt times
        un = u.copy()  # copy the existing values of u into un
        iterations.append(un)
        for i in range(1, len(cb)-2):
            u[i] = un[i] + dt * (((nu / (dx[i] ** 2)) * ((un[i] ** 2) + 2 * un[i] * un[i-1] + (un[i-1] ** 2))) - ((c / dx[i]) * (un[i] - un[i-1])))
    with open('output/%s.csv' %filename, 'w', newline='') as f:
        writer = csv.writer(f)
        writer.writerows(iterations)

Klar müsste ich da einiges ändern um das zum Laufen zu bekommen aber man nehme mal an ich möchte nur den numerischen Teil überführen mit der Berechnung. Wie müsste ich da vorgehen? Ist das überhaupt eine gute Idee? Gibt es dafür bessere Lösungen?

Momentan stecke ich hier ziemlich am Beginn und weiß schlicht nicht wie ich das ansetzen soll bzw. wie ich dann damit rechnen kann/soll.

class Cell:

    def __init__(self, u):
        self.u = 0.0
       
c1 = Cell('0')
c2 = Cell('1')
c3 = Cell('0')

Am liebsten wäre mir sowas wie: Variable u neu der Zelle n = Variable U der Zelle n + Zelle n+1 mit Variable u

Danke schonmal und liebe Grüße!

P.s. Bevor jemand fragt, nein das sind keine Hausaufgaben. Ich spiel schon länger mit Python rum aber habe mir bisher nicht an die Objektorientierung rangetraut da ich diese im Studium auch schon nicht verstanden habe..^^

Ich hätte auch noch dieses Codeschnipsel hier gefunden das ja theoretisch im Kontext meines Hauptprogramms funktioniert, allerdings checke ich nicht ganz was es macht. Den Ersteller kann ich leider nicht mehr fragen...

cells = []

for j in range(len(u)):
    cell = Cell()   # Was zum Geier macht der Teil des Codes?
    cell._index = j
    cell._midpoint = cm[j]
    if len(cells) > 0:
        cells[-1]._right_ngb = cell

    cells.append(cell)
#     print(cm[j])

len(u) wäre die Anzahl an Elementen im Vektor u.
cm = Zellenmittlepunkt
right ngb = right neighbour

 

[DeusExMachina]

Doofe frage warum eine Klasse?

cell = Cell() # Was zum Geier macht der Teil des Codes?

Erzeugt eine Instanz der Class Cell und weißt sie der Variablen cell zu.


iterations = [] # neu laden der Variable um alte Werte loszuwerden

Ähm nein. Nicht neu laden sondern mit einem leeren Array überschreiben.

 

[DeusExMachina]

Ich hab noch mal auf dem rum gedacht was du da versuchst zu bauen.
Und so wirklich verstehe ich es nicht.


Wenn du ein Gitter aufspannen willst solltest du dir erstmal über die Struktur im klaren sein.
Ganz Simple könntest du einen 2D Array nehmen
wobei der erste Index die Spalte der zweite die Zeile darstellt

also so was gridColumns = [ [A, B, C], [A, B, C], [A, B, C] ] wobei dann A, B, C jeweils Instanzen deiner Cell Class sind.

Aber wie gesagt ich verstehe nicht ganz was du machen willst

 

[Tenferenzu]

Ganz Simple könntest du einen 2D Array nehmen

Ist 1D nicht einfach genug?^^

Ich hab noch mal auf dem rum gedacht was du da versuchst zu bauen.
Und so wirklich verstehe ich es nicht.

Ich möchte zu Beginn mal einen Impuls (also eine Eins neben vielen Nullen) numerisch durch das Grid wandern lassen. Mit dem aktuellen Code klappt das auch ganz gut, allerdings sind halt eine ziemlich große Anzahl an Schleifen involviert. Von der Klasse würde ich mir erhoffen, das bissl einschränken und übersichtlicher machen zu können.

Das hier wäre der Output vom aktuellen Code wenn ich einen Entsprechenden Impuls mit einer Breite von 1 und einer Höhe von 1 zu Beginn aufgebe.

Verwendet wird hierfür die ganz allgemeine Transportgleichung in numerischer Form mit einem angepassten Upwind Schema.

Ich hoffe das Beantwortet mehr Fragen als es jetzt ev. noch neu dazugekommen sind

 

[DeusExMachina]

Ein normales eindimensionales Array würde reichen wenn du kein Grind haben wolltest.

zum Beispiel wenn du dich entlang der X-Achse bewegen willst und jeder Eintrag einen X-Wert repräsentiert.
Nehmen wir hier mal als Beispiel die 20.
Dich interessiert dann ja nur der Y-Wert zum X

Also könnte deine Klasse zum Beispiel so aussehen:

class Point:
        
    def __init__(self, x):
        self.x = x
        self.y = self.calc_y(x)
    
    def calc_y(x):
        #do some calculations
        return y

wenn du dann ein neues Object für eine Celle mit cell = Cell(x) erstellst solltest du mit cell.y den wert bekommen

 

[simpsonsfan]

Kann gerade nicht auf alles eingehen, später evtl. detaillierter.
Zunächst mal ist es performancetechnisch eher keine gute Idee, die Zellen eines FV-Lösers als High-Level Objekte zu betrachten. Dort will man möglichst vektorisierte Befehle (SSE/AVX) verwenden. Zu Übungszwecken aber natürlich vollkommen i.O.

Schleifen wirst du auch brauchen, wenn deine Zellen eine eigene Klasse haben. Nur iterierst du dann halt über alle deine Zell-Instanzen. Ggf. kannst du den Code mit der Objektorientierung übersichtlicher gestalten. Vielleicht suchst du dir auch besser ein anderes Beispiel zu Beginn.

 

[Tenferenzu]

Zunächst mal ist es performancetechnisch eher keine gute Idee, die Zellen eines FV-Lösers als High-Level Objekte zu betrachten. Dort will man möglichst vektorisierte Befehle (SSE/AVX) verwenden.

Genau das ist eine meiner Befürchtungen. Ich bin jetzt zwar kein Experte aber mit numpy ist das alles schon um eine Größenordnung schneller als ohne. Bei 10 Zellen und ein paar Zeitschritten macht das nicht viel aus aber bei 1-100k Zellen und 2-10k Zeitschritten ist der PC schon ordentlich am Rechnen.

Zu Übungszwecken aber natürlich vollkommen i.O.

Danke!

Ggf. kannst du den Code mit der Objektorientierung übersichtlicher gestalten.

Das wäre eines der Ziele. Ich weiß nur nicht ganz ob das von Grund auf so ansetzen sollte oder erst sobald man alle Routinen etc. mal fertig hat. Aktuell tendiere ich ja zu letzterem.

Vielleicht suchst du dir auch besser ein anderes Beispiel zu Beginn.

Hast du Empfehlungen? Gerne auch mit etwas komplexeren Rechenbeispielen. Die meisten Tutorials arbeiten meist nur mit Personen, Namen, etc. und geben dann einen Satz oder ähnliches aus.

Kann gerade nicht auf alles eingehen, später evtl. detaillierter.

Danke schonmal!

@DeusExMachina Danke! Ich schau mir das mal an wenn ich später zu Hause bin.

 

[DeusExMachina]

Kein Ding.
Das ist aber auch keine echte Lösung und Performance habe ich mit Absicht überhaupt nicht berücksichtigt.

Mir ging es erstmal nur um die Abbildung eines Punktes als Objekt

 

[simpsonsfan]

So, noch ein bisschen mehr dazu.
Erstmal die Frage, möchtest du nur eine skalare Größe (in deine Plot und deiner Gleichungsangabe bspw. eine Konzentration c) durch deine Domain mit der fixen konvektiven Geschwindigkeit u transportieren? Oder möchtest du einen Impuls transportieren (und damit eine Geschwindigkeit)
Wenn ich das recht sehe, dann hast du in Zeile 20 einen Fehler, nämlich u und c vertauscht.
Bei dir ist dann u die Transportgröße und c die konstante konvektive Geschwindigkeit.
Es wäre mMn schöner, man würde hier das c durch ein φ ersetzen. Das ist klarere Konvention für ein beliebiges Skalar, während c nach belieben für eine Konzentration, eine Geschwindigkeit oder beliebige Konstanten stehen kann.
Und da bei Strömungen eben auch der Impuls transportiert wird, ist die räumlich und zeitlich variable Geschwindigkeit oft gesucht, so dass deine Transportgröße die Geschwindigkeit (Impuls/Masse) u ist. Daher wird in solchen Codes oft einfach u als Symbol verwendet, auch wenn man eine andere Größe transportieren möchte.
Ich sehe da zumindest Potential zur Verwirrung.
Das "n-1" in deiner Gleichungsangabe (bei der zweiten Ableitung) müsste dann auch wieder ein Tippfehler sein und nur "n" sein, oder?

Ich weiß nur nicht ganz ob das von Grund auf so ansetzen sollte oder erst sobald man alle Routinen etc. mal fertig hat. Aktuell tendiere ich ja zu letzterem.

Ich würde eher zum ersten tendieren. Die gesamte Struktur deines Codes wird eine andere. (Und zwar insgesamt eher eine übersichtlichere.)
Du solltest dir am besten vorher überlegen, wie eine Kapselung Sinn macht.
Hier hast du bspw. deine Zellen, die können eine Klasse werden.
Jede Zelle hat dann eine Eigenschaft "skalarer Wert φ." Eventuell auch die Information "linker Nachbar" und "rechter Nachbar." Ggf. noch "Zellweite."
Jede Zelle kann eine Methode "gehe einen Timestep vorwärts/update()" ausführen.
Die Methode "update()" kann später möglicherweise verschiedene Rechenschemas verwenden, bspw. FirstOrderUpwind oder SecondOrderUpwind.

Dann hast du ein Gitter. Das wird durch eine Klasse abgebildet und enthält Verweise auf alle Zellen.
Das Gitter (oder weiter getrennt/abstrahiert ein Solver) hat bspw. die Methoden "initialisiereFeld(Dateiname)" und "löseFeld(timestep, endtime, velocity, diffusivity)."
"löseFeld" wiederum ruft dann iterativ (in jedem Timestep) die update()-Methode jeder Zelle auf. Also bspw.

def löseFeld:
  time t=0
  while (t< endtime):
    solver.step()
    t += dt

def step:
  for cell in gitter.cells:
    cell.update()
    solver.solutionhistory.appendCurrentIteration()

Für die Datenstruktur des zeitlichen Verlaufs musst du dir überlegen, wie das gespeichert wird. Erhält jede Zelle einen zeitlichen Verlauf ihrer Lösung? Oder erhält doch besser die Solver-Klasse neben dem Gitter noch eine Feldlösung zu beliebigen Zeiten? Dann kannst du später (nachdem das Programm gelaufen ist) über solver.getSolution(Timestep = 37) die Lösung zum Zeitpunkt 37 abrufen.

Das Schöne ist, dass du damit den Code besser strukturierst und Änderungen dann später hoffentlich einigermaßen zügig einbauen kannst, weil du davor eben schonmal gekapselt hast und überlegt hast, welche Daten wohin müssen und welche Objekte auf was zugreifen müssen.

Ist jetzt etwas viel Text, aber vielleicht hilft es dir ja ein bisschen.

 

[Tenferenzu]

Danke für deinen Input @simpsonsfan, im Endeffekt wird es wohl darauf hinauslaufen, dass ich den SIMPLE Algorithmus implementieren werde um die Druck-, Geschwindigkeitskopplung hinzubekommen. So viel Auswahl an Algorithmen gibt es ja nicht um das zu realisieren was ich machen will.
Danke nochmal! Du hast mir echt weitergeholfen!

 

[simpsonsfan]

Gerne.
Zum Thema SIMPLE: Was genau möchtest du denn machen? Es gibt ja durchaus einige verschiedene Methoden zur numerischen Lösung der Navier-Stokes-Gleichungen (bzw. Euler-Gleichungen.)
Den segregated-Ansatz habe ich in den Vorlesungen eigentlich immer als etwas weniger intuitiv empfunden, als eben bspw. gleich gekoppelt zu rechnen und eben ein Riemannproblem an jeder Face/Zellkante zu lösen. Das dann explizit, dann kann man einigermaßen nachvollziehen, was passiert.
Es muss also kein SIMPLE oder überhaupt ein Druckkorrekturverfahren sein. (Sofern das nicht explizit dein Ziel ist; Und natürlich kann es das auch ruhig sein - ich wollte nur auf die Aussage antworten, dass es nicht so viel Auswahl gibt. Ein bisschen zumindest schon.)

 

[Tenferenzu]

Was genau möchtest du denn machen?

Ganz genau? Ich möchte elektrochemische Membranprozesse simulieren können. Dazu brauche ich zumindest mal den Elektrochemieteil, das Lösungsdiffusionsmodell und das Porenflussmodell für den Teil der Membran der für die mechanische Stabilität verantwortlich ist.

Innerhalb der Membran kann ich mir relativ gut mit der Diffusion/dem chemischen Potential abhelfen. Die Schwierigekeit für mich persönlich ist dann der Fluss im Lumen (dem inneren der Hohlfasermembran bzw. der Feedseite) und dem Fluss auf der Permeatseite.

Es muss also kein SIMPLE oder überhaupt ein Druckkorrekturverfahren sein. (Sofern das nicht explizit dein Ziel ist; Und natürlich kann es das auch ruhig sein - ich wollte nur auf die Aussage antworten, dass es nicht so viel Auswahl gibt. Ein bisschen zumindest schon.)

Mir würde es in erster Linie um die Grenzschichten um die Membran als auch die Verteilung innerhalb der Membran gehen. Der Rest drumherum ist nur Mittel zum Zweck um eine Spezies dorthin transportieren zu können bzw. auf der Permatseite diese wieder abtransportieren zu können. Wenn ich dort nun sowieso annehme, dass keine Turbulenz herrscht, hätte ich horizontal (laut Skizze) sowieso nur Diffusion. Wo ich dann wirklich zu kämpfen habe, ist der Abtransport in die Membran. Klar, direkt an der Membran ist die Konzentration von XYZ jetzt nicht so hoch aber es müsste ja auch noch was nachfließen da ein Teil des Gases in der Membran gelöst wird -> Druchverlust und zumindest etwas Strömung in richtung Membran.

Klar gibt's dann noch den Druckverlust innerhalb der Hohlfaser da diese ja auch entsprechend lang und dünn sein kann. Der Druck geht dann wieder in die Membranmodelle ein.

Hier mal eine kleine Skizze wie ich mir das so ca. vorstelle.

Wenn das ganze viel einfacher zu lösen wäre, würde ich mich nicht beschweren wenn du da ev. Tipps hast

Ich bin SOOO knapp davor das als UDF in Fluent umzusetzen aber der innere Schweinehund will es FOSS halten und lieber was haben das überall funktioniert.

Vorlesungen

Das Problem ist, die fehlen mir halt^^ Ich muss mir das gerade selber beibringen. Wir haben zwar mit Fluent und OpenFOAM gearbeitet aber das waren Tutorials/Aufgaben die 'relativ' leicht waren. Numerik in dem Sinne gab es nicht in meinem Studienplan deshalb ist das gerade etwas viel

 

[simpsonsfan]

Fluent und seine UDFs haben auch irgendwo ein gewöhnungsbedürftiges Handling.
Wenn FOSS das Kriterium ist, würde ich aber eher bei OpenFOAM bleiben. Das lässt sich eigentlich ganz gut anpassen. Würde mich auch nicht wundern, falls es da irgendwo schon Ansätze in die Richtung gibt.

Grundsätzlich könnte es wahrscheinlich auch nicht schaden, sich zumindest schonmal bei den kommerziellen Anbietern umzuschauen. Richtung Brennstoffzellen und Batterien tut sich bei StarCCM in den letzten Releases auch bisschen was. Hab das aber selber noch nicht wirklich angeschaut.
Aber bevor du da was in Python/Numpy wie auch immer from Scratch aufsetzt, würde ich lieber OpenFOAM nehmen, mit der Option, das dann später ggf. auch 3D mit Originalgeometrie rechnen zu können.