From f6150c7682ad2d77e124c5f6eab1d6f069f437e5 Mon Sep 17 00:00:00 2001
From: Jussi Aittoniemi <jussi.aittoniemi@bluware.com>
Date: Tue, 7 Nov 2023 14:26:56 +0100
Subject: [PATCH] Fix crust properties per-node

---
 subsheat3D/Helpers.py |  1 +
 warmth/mesh_model.py  | 33 ++++++++++++++++-----------------
 2 files changed, 17 insertions(+), 17 deletions(-)

diff --git a/subsheat3D/Helpers.py b/subsheat3D/Helpers.py
index 7ce037b..2b76823 100644
--- a/subsheat3D/Helpers.py
+++ b/subsheat3D/Helpers.py
@@ -25,6 +25,7 @@ class NodeParameters1D:
     kCrust: float = 2.5
     kAsth: float = 100
     rhp: float = 2
+    crustRHP: float = 2
     crustliquid: float = 2500.0
     crustsolid: float = 2800.0
     lithliquid: float = 2700.0
diff --git a/warmth/mesh_model.py b/warmth/mesh_model.py
index ef524b1..b0d76d4 100644
--- a/warmth/mesh_model.py
+++ b/warmth/mesh_model.py
@@ -9,7 +9,7 @@
 from warmth.build import single_node
 from .model import Model
 from warmth.logging import logger
-from .mesh_utils import  NodeParameters1D, top_crust,top_sed,thick_crust,  top_lith, top_asth, top_sed_id, bottom_sed_id,NodeGrid
+from .mesh_utils import  top_crust,top_sed,thick_crust,  top_lith, top_asth, top_sed_id, bottom_sed_id,NodeGrid
 from .resqpy_helpers import write_tetra_grid_with_properties, write_hexa_grid_with_properties,read_mesh_resqml_hexa
 def tic():
     #Homemade version of matlab tic and toc functions
@@ -79,11 +79,8 @@ def __init__(self, model:Model, modelName="test", sedimentsOnly = False):
         self.mean_porosity = None
         self.c_rho = None
         self.numberOfSediments = model.builder.input_horizons.shape[0]-1 #skip basement
-        # self.globalSediments = self.node1D[0].sediments.copy()
 
-        self.parameters1D = self.node1D[0].parameters if hasattr(self.node1D[0], 'parameters') else NodeParameters1D()
         self.interpolators = {}
-        print("Using 1D Node parameters", self.parameters1D)
    
     def write_tetra_mesh_resqml( self, out_path):
         """Prepares arrays and calls the RESQML output helper function:  the lith and aesth are removed, and the remaining
@@ -315,17 +312,17 @@ def cRhoForLayerID(self, ss, node_index):
         #
         # prefactor 1000 is the heat capacity.. assumed constant
         #
+        node = self.node1D[node_index]
         if (ss==-1):
-            return 1000*self.parameters1D.crustsolid
+            return 1000*node.crustsolid
         if (ss==-2):
-            return 1000*self.parameters1D.lithsolid
+            return 1000*node.lithsolid
         if (ss==-3):
-            return 1000*self.parameters1D.crustsolid
+            return 1000*node.crustsolid
         if (ss>=0) and (ss<self.numberOfSediments):
-            node = self.node1D[node_index]
             rho = node.sediments.solidus[ss]
             return 1000*rho
-        return 1000*self.parameters1D.crustsolid
+        return 1000*node.crustsolid
 
     def kForLayerID(self, ss, node_index):
         """Thermal conductivity for a layer ID index
@@ -339,19 +336,19 @@ def kForLayerID(self, ss, node_index):
         #     print(len(self.cell_data_layerID), np.amin(self.cell_data_layerID), np.amax(self.cell_data_layerID))
         #     print("kForLayerID", ss, ind0, ind1, ind2, self.cell_data_layerID[self.cell_index[cell_id]] )
         # assert self.cell_data_layerID[ind2] == ss
+        if (node_index > len(self.node1D)-1):
+            print("cell ID", node_index, len(self.node1D))
+            breakpoint()
+        node = self.node1D[node_index]
         if (ss==-1):
-            return self.parameters1D.kCrust
+            return node.kCrust
         elif  (ss==-2):
-            return self.parameters1D.kLith
+            return node.kLith
         elif (ss==-3):
-            return self.parameters1D.kAsth
+            return node.kAsth
         elif (ss>=0) and (ss<self.numberOfSediments):
             # kc = self.globalSediments.k_cond[ss]
             # node_index = cell_id // (self.numberOfSediments+6)  # +6 because crust, lith, aest are each cut into two
-            if (node_index > len(self.node1D)-1):
-                print("cell ID", node_index, len(self.node1D))
-                breakpoint()
-            node = self.node1D[node_index]
             kc = node.sediments.k_cond[ss]
             return kc
 
@@ -360,7 +357,9 @@ def rhpForLayerID(self, ss, node_index):
         """Radiogenic heat production for a layer ID index
         """
         if (ss==-1):
-            return 0
+            node = self.node1D[node_index]
+            kc = node.crustRHP
+            return kc
         elif (ss==-2):
             return 0
         elif (ss==-3):