Input: Add TopoHeight to the interface (#123)

* Input: Add TopoHeight to the interface * add conditional when there is no free_surf * deal with the memory, free topo_height * add some conditionals to prevent segfault * valgrind fixes * use topo instead of topo_chain * add new Setup * fix Main_DOODZ * need more input from Roman... * parse integrated boundary flux instead of altitude * Update Setup.c * No major change, just relabeling * correct value * add new zero mean topo function and rename structs * rename setup file * activate zero mean topo in RiverTomLowerCrust * force adding headers * add KeepZeroMeanTopo to header * ready * this should fix RiverTom * Ready to test? --------- Co-authored-by: tduretz <[email protected]> Co-authored-by: Tgeffroy <[email protected]> Co-authored-by: Thibault Duretz <[email protected]> Co-authored-by: thibault duretz <[email protected]>
tduretz · Mar 11, 2024 · 08544b3 · 08544b3
1 parent fa60e86
commit 08544b3
Show file tree

Hide file tree

Showing 9 changed files with 440 additions and 19 deletions.
diff --git a/MDLIB/FreeSurface.c b/MDLIB/FreeSurface.c
@@ -412,6 +412,7 @@ void AllocateMarkerChain( surface *topo, markers* topo_chain, params model ) {
     topo_chain->Vx          = DoodzCalloc( topo_chain->Nb_part_max, sizeof(DoodzFP) );
     topo_chain->Vz          = DoodzCalloc( topo_chain->Nb_part_max, sizeof(DoodzFP) );
     topo_chain->phase       = DoodzCalloc( topo_chain->Nb_part_max, sizeof(int) );
+    // for (int i=1; i<topo_chain->Nb_part_max; i++) topo_chain->phase[i] = -1; 
     topo->height            = DoodzCalloc( (model.Nx),sizeof(DoodzFP) );
     topo->height0           = DoodzCalloc( (model.Nx),sizeof(DoodzFP) );
     topo->vx                = DoodzCalloc( (model.Nx),sizeof(DoodzFP) );
@@ -1051,3 +1052,27 @@ void SurfaceVelocity( grid *mesh, params model, surface *topo, markers* topo_cha
 /*--------------------------------------------------------------------------------------------------------------------*/
 /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
 /*--------------------------------------------------------------------------------------------------------------------*/
+
+void KeepZeroMeanTopo(params *model, surface *topo, markers *topo_chain ) {
+    printf("Make zero mean topo\n");
+    double mean_z = 0.0;
+    for (int i=0; i<model->Nx; i++){
+        mean_z += topo->height[i];
+    }
+    mean_z /= model->Nx;
+    for (int i=0; i<model->Nx; i++) {
+        topo->height[i] -= mean_z;
+    }
+    mean_z = 0.0;
+    int nmark = 0;
+    for (int i=0; i<topo_chain->Nb_part_max; i++){
+        if (topo_chain->phase[i] > -1) {
+        mean_z += topo_chain->z[i];
+        nmark  ++;
+        }
+    }
+    mean_z /= nmark;
+    for (int i=0; i<topo_chain->Nb_part_max; i++){
+        topo_chain->z[i] -= mean_z;
+    }
+}
diff --git a/MDLIB/InputOutput.c b/MDLIB/InputOutput.c
@@ -1062,6 +1062,7 @@ Input ReadInputFile( char *fileName ) {
     model.zmin               = ReadDou2( fin, "zmin",    -1.0 )/scaling.L;  // Spatial domain extent  
     model.zmax               = ReadDou2( fin, "zmax",     1.0 )/scaling.L;  // Spatial domain extent
     model.balance_boundaries = ReadInt2( fin, "balance_boundaries",   0 );  // Switch this to activate boundary velocity balancing
+    model.zero_mean_topo     = ReadInt2( fin, "zero_mean_topo",       0 );  // Switch this to force zero mean topography
     // Time domain
     model.Nt                 = ReadInt2( fin, "Nt",         1 );            // Number of time steps    
     model.dt                 = ReadDou2( fin, "dt",       0.0 ) /scaling.t; // Time step

diff --git a/MDLIB/Main_DOODZ.c b/MDLIB/Main_DOODZ.c
@@ -52,13 +52,21 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
     asprintf(&BaseOutputFileName, "Output");
     asprintf(&BaseParticleFileName, "Particles");
 
-    MdoodzInput input = (MdoodzInput) {
-      .model = inputFile.model,
-      .scaling = inputFile.scaling,
-      .materials = inputFile.materials,
-      .crazyConductivity = NULL,
+    MdoodzInput input = (MdoodzInput){
+            .model             = inputFile.model,
+            .scaling           = inputFile.scaling,
+            .materials         = inputFile.materials,
+            .crazyConductivity = NULL,
+            .flux       = NULL,
     };
 
+    if (input.model.free_surface) {
+      input.flux = malloc(sizeof(LateralFlux));
+      if (input.flux != NULL) {
+        *input.flux = (LateralFlux){.east = -0.0e3, .west = -0.0e3};
+      }
+    }
+
     if (setup->MutateInput) {
       setup->MutateInput(&input, setup->mutateInputParams);
     }
@@ -103,7 +111,7 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
     markers particles = PartAlloc(inputFile.particles, &input.model );
 
     // Allocate marker chain
-    markers       topo_chain, topo_chain_ini;
+    markers      topo_chain, topo_chain_ini;
     surface      topo, topo_ini;
     if (input.model.free_surface == 1 ) AllocateMarkerChain( &topo,     &topo_chain, input.model );
     if (input.model.free_surface == 1 ) AllocateMarkerChain( &topo_ini, &topo_chain_ini, input.model );
@@ -128,7 +136,8 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
         }
 
         // Initial grid tags
-        SetBCs(*setup->SetBCs, &input, &mesh);
+        SetBCs(*setup->SetBCs, &input, &mesh, &topo);
+
         if (input.model.free_surface == 1 ) {
 
             // Define the horizontal position of the surface marker chain
@@ -190,7 +199,7 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
           printf("Running with crazy conductivity for the asthenosphere!!\n");
         }
         // Initial solution fields
-        SetBCs(*setup->SetBCs, &input, &mesh);
+        SetBCs(*setup->SetBCs, &input, &mesh, &topo);
 
         if (input.model.mechanical == 1) {
             InitialiseSolutionFields( &mesh, &input.model );
@@ -216,7 +225,7 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
         P2Mastah( &input.model, particles, input.materials.Cv,    &mesh, mesh.Cv, mesh.BCp.type,  0, 0, interp, cent, 1);
         P2Mastah( &input.model, particles, input.materials.Qr,    &mesh, mesh.Qr, mesh.BCp.type,  0, 0, interp, cent, 1);
 
-        SetBCs(*setup->SetBCs, &input, &mesh);
+        SetBCs(*setup->SetBCs, &input, &mesh, &topo);
         if (input.model.initial_cooling == 1 ) ThermalSteps( &mesh, input.model,  mesh.rhs_t, &particles, input.model.cooling_duration, input.scaling );
         if (input.model.therm_perturb == 1 ) SetThermalPert( &mesh, input.model, input.scaling );
         Interp_Grid2P_centroids2( particles, particles.T,    &mesh, mesh.T, mesh.xvz_coord,  mesh.zvx_coord,  mesh.Nx-1, mesh.Nz-1, mesh.BCt.type, &input.model );
@@ -584,7 +593,7 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
                 printf("Running with normal conductivity for the asthenosphere...\n");
             }
             // Allocate and initialise solution and RHS vectors
-            SetBCs(*setup->SetBCs, &input, &mesh);
+            SetBCs(*setup->SetBCs, &input, &mesh, &topo);
 
             // Reset fields and BC values if needed
             //        if ( input.model.pure_shear_ALE == 1 ) InitialiseSolutionFields( &mesh, &input.model );
@@ -1019,6 +1028,8 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
                     CorrectTopoIni( &particles, input.materials, &topo_chain_ini, &topo, input.model, input.scaling, &mesh);
                     MarkerChainPolyFit( &topo_ini, &topo_chain_ini, input.model, mesh );
 
+                    if ( input.model.zero_mean_topo == 1 ) KeepZeroMeanTopo( &input.model, &topo, &topo_chain );
+
                     // Sedimentation
                     if ( input.model.surface_processes == 2 ) {
                         AddPartSed( &particles, input.materials, &topo_chain, &topo, input.model, input.scaling, &mesh);
@@ -1103,6 +1114,16 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
             DoodzFree( Jacob.eqn_p );
         }
 
+        // double mean_z = 0.;
+        // int n=0;
+        // for (int i=0; i<topo_chain.Nb_part_max; i++) {
+        //     if (topo_chain.phase[i] != -1) {
+        //         mean_z += topo_chain.z[i];
+        //         n++;
+        //     }
+        // }
+        // printf("mean topo on markers = %2.2e --- %d out of %d and max is %d \n", mean_z/n*input.scaling.L, n, topo_chain.Nb_part, topo_chain.Nb_part_max);
+
         //------------------------------------------------------------------------------------------------------------------------------//
 
         // Write output data
@@ -1157,6 +1178,10 @@ void RunMDOODZ(char *inputFileName, MdoodzSetup *setup) {
     // Free arrays
     GridFree( &mesh, &input.model );
 
+    if (input.flux) {
+      free(input.flux);
+    }
+
     if (input.crazyConductivity) {
         free(input.crazyConductivity->phases);
         free(input.crazyConductivity);

diff --git a/MDLIB/RheologyDensity.c b/MDLIB/RheologyDensity.c
@@ -452,7 +452,7 @@ double ViscosityConcise( int phase, double G, double T, double P, double d, doub
   *d1 = materials->gs_ref[phase];
 
   // Tensional cut-off
-  if ( model->gz<0.0 && P<0.0     ) { P = 0.0; } // printf("Aie aie aie P < 0 !!!\n"); exit(122);
+  //if ( model->gz<0.0 && P<0.0     ) { P = 0.0; } // printf("Aie aie aie P < 0 !!!\n"); exit(122);
 
   // Visco-plastic limit
   if ( elastic==0                 ) { G = 1e1; dil = 0.0; K = 1e10;}; //K = 1e1;
@@ -689,7 +689,7 @@ double ViscosityConcise( int phase, double G, double T, double P, double d, doub
 
       // In case return mapping has failed (because of tension), return to a von Mises minimum stress 
       if (Tiic<0.0) {
-        printf("Aie, tension!\n"); exit(1);
+        //printf("Aie, tension!\n"); exit(1);
         F_trial = Tii - Tiimin;
         gdot    = F_trial /  (eta_ve);
         Tiic    = Tii - eta_ve*gdot;

diff --git a/MDLIB/Setup.c b/MDLIB/Setup.c
@@ -176,7 +176,7 @@ void ValidateInternalPoint(POSITION position, char bcType, Coordinates coordinat
 /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
 /*--------------------------------------------------------------------------------------------------------------------*/
 
-void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh) {
+void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh, surface *topo) {
   /* --------------------------------------------------------------------------------------------------------*/
   /* Set the BCs for Vx on all grid levels */
   /* Type  0: Dirichlet point that matches the physical boundary (Vx:
@@ -195,7 +195,6 @@ void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh) {
   double VxWestSum = 0.0;
   double VxEastSum = 0.0;
 
-
   for (int l = 0; l < mesh->Nz + 1; l++) {
     for (int k = 0; k < mesh->Nx; k++) {
       const int c = k + l * (mesh->Nx);
@@ -387,6 +386,12 @@ void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh) {
     free(newBoundary);
   }
 
+  // override by integrated velocity values
+  if (instance->flux != NULL) {
+    instance->flux->west = VxWestSum*mesh->dz;
+    instance->flux->east = VxEastSum*mesh->dz;
+  }
+
   /* --------------------------------------------------------------------------------------------------------*/
   /* Set the BCs for Vz on all grid levels */
   /* Type  0: Dirichlet point that matches the physical boundary (Vx:
@@ -458,8 +463,8 @@ void SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh) {
     }
   }
 
-  printf("VzWestSum: %f, VzEastSum: %f: \n", VzWestSum, VzEastSum);
-  printf("total West+East+South sum: %f\n", VxWestSum + VxEastSum - VzSouthSum);
+  printf("VxWestSum*dx: %f, VxEastSum*dx: %f, VzEastSum*dz: %f: \n", VxWestSum*mesh->dz, VxEastSum*mesh->dz, VzSouthSum*mesh->dx);
+  printf("Total West+East+South sum: %f\n", fabs(VxWestSum*mesh->dz) + fabs(VxEastSum*mesh->dz)  - fabs(VzSouthSum*mesh->dx) );
 
   if (instance->model.balance_boundaries && (VzSouthSum > tolerance || VzSouthSum < -tolerance)) {
     int zeroValuesCount = 0;
@@ -994,4 +999,4 @@ bool IsRectangleCoordinates(Coordinates coordinates, Rectangle rectangle, double
 
 /*--------------------------------------------------------------------------------------------------------------------*/
 /*------------------------------------------------------ M-Doodz -----------------------------------------------------*/
-/*--------------------------------------------------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------------------------------------------------*/
diff --git a/MDLIB/include/mdoodz.h b/MDLIB/include/mdoodz.h
@@ -42,7 +42,7 @@ typedef struct {
 
 // params contains the model parameters
 typedef struct {
-  int    balance_boundaries;
+  int    balance_boundaries, zero_mean_topo; 
   char   description[500];
   double xmin, zmin, xmax, zmax, time, dx, dz, dt, dt0, dt_start, dt_max, L0,
           dt_min;
@@ -297,11 +297,17 @@ struct MdoodzSetup {
   MutateInputParams         *mutateInputParams;
 };
 
+typedef struct {
+  double east;
+  double west;
+} LateralFlux;
+
 struct MdoodzInput {
   char              *inputFileName;
   params             model;
   mat_prop           materials;
   scale              scaling;
+  LateralFlux        *flux;
   CrazyConductivity *crazyConductivity;
   Geometry          *geometry;
 };

diff --git a/MDLIB/mdoodz-private.h b/MDLIB/mdoodz-private.h
@@ -1,6 +1,7 @@
 #include "cholmod.h"
 #include "cs.h"
 #include "mdoodz.h"
+
 //---------------------------------------------------------------------------------------------------------------------------//
 //------------------------------------------------ MACROS DEFINITIONS -------------------------------------------------------//
 //---------------------------------------------------------------------------------------------------------------------------//
@@ -419,6 +420,7 @@ void            UpdateDensity(grid *, markers *, mat_prop *, params *, scale *);
 void            DiffuseAlongTopography(grid *, params, scale, double *, double *, int, double, double);
 void            AddPartSed(markers *, mat_prop, markers *, surface *, params, scale, grid *);
 void            CorrectTopoIni(markers *, mat_prop, markers *, surface *, params, scale, grid *);
+void            KeepZeroMeanTopo(params*, surface*, markers* );
 
 // Decoupled solver
 void            KillerSolver(SparseMat *, SparseMat *, SparseMat *, SparseMat *, DirectSolver *, double *, double *, double *, params, grid *, scale, SparseMat *, SparseMat *, SparseMat *, SparseMat *, SparseMat *);
@@ -573,7 +575,7 @@ void            MinMaxArray(double *array, double scale, int size, char *text);
 
 void            BuildInitialTopography(BuildInitialTopography_ff buildInitialTopography, MdoodzInput *instance, markers *topo_chain);
 void            SetParticles(SetParticles_ff setParticles, MdoodzInput *instance, markers *particles);
-void            SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh);
+void            SetBCs(SetBCs_ff setBCs, MdoodzInput *instance, grid *mesh, surface *topo);
 
 void            ValidateSetup(MdoodzSetup *setup, MdoodzInput *instance);