parameters.hpp

#ifndef DYNEARTHSOL3D_PARAMETERS_HPP
#define DYNEARTHSOL3D_PARAMETERS_HPP

#include <map>
#include <string>
#include <utility>
#include <vector>
#include <unordered_map>

#include <nanoflann.hpp>
#ifdef NPROF
#include <nvtx3/nvToolsExt.h> 
// #include <nvToolsExt.h> 
#endif

#include "constants.hpp"
#include "array2d.hpp"
#ifdef ACC
#include "knn_bvh.hpp"
#endif

// Forward declarations
#ifdef HAS_GOSPL_CPP_INTERFACE
class GoSPLDriver;
#endif

typedef std::pair<int,int> int_pair;
typedef std::pair<double,double> double_pair;
typedef std::unordered_map<int,int> int_map;
typedef std::vector<int_map> int_map2D;

typedef std::vector<double> double_vec;
typedef std::vector<int> int_vec;
typedef std::vector<int_vec> int_vec2D;
typedef std::vector<double_vec> double_vec2D;
typedef std::vector<uint> uint_vec;
typedef std::vector<unsigned char> uchar_vec;
typedef std::vector<bool> bool_vec;
typedef std::vector<size_t> size_t_vec;
typedef std::vector<int_pair> int_pair_vec;

typedef Array2D<double,NDIMS> array_t;
typedef Array2D<double,NSTR> tensor_t;
typedef Array2D<double,NODES_PER_ELEM> shapefn;
typedef Array2D<double,1> regattr_t;
typedef Array2D<double,NODES_PER_ELEM*3> elem_cache;

typedef Array2D<int,NODES_PER_ELEM> conn_t;
typedef Array2D<int,NDIMS> segment_t;
typedef Array2D<int,1> segflag_t;
typedef Array2D<int,NODES_PER_CELL> regular_t;
typedef nanoflann::KNNResultSet<double> KNNResultSet;

typedef array_t::Accessor ArrayAccessor;
typedef tensor_t::Accessor TensorAccessor;
typedef shapefn::Accessor ShapefnAccessor;
// typedef regattr_t::Accessor RegattrAccessor;
typedef elem_cache::Accessor ElemCacheAccessor;

typedef conn_t::Accessor ConnAccessor;
typedef segment_t::Accessor SegmentAccessor;
// typedef segflag_t::Accessor SegflagAccessor;
// typedef regular_t::Accessor RegularAccessor;

typedef array_t::ConstAccessor ConstArrayAccessor;
typedef tensor_t::ConstAccessor ConstTensorAccessor;
typedef shapefn::ConstAccessor ConstShapefnAccessor;
// typedef regattr_t::ConstAccessor ConstRegattrAccessor;
typedef elem_cache::ConstAccessor ConstElemCacheAccessor;

typedef conn_t::ConstAccessor ConstConnAccessor;
typedef segment_t::ConstAccessor ConstSegmentAccessor;
// typedef segflag_t::ConstAccessor ConstSegflagAccessor;
typedef regular_t::ConstAccessor ConstRegularAccessor;

typedef array_t::ConstIndirectAccessor ConstArrayIndirectAccessor;
// typedef tensor_t::ConstIndirectAccessor ConstTensorIndirectAccessor;
// typedef shapefn::ConstIndirectAccessor ConstShapefnIndirectAccessor;
// typedef regattr_t::ConstIndirectAccessor ConstRegattrIndirectAccessor;
// typedef elem_cache::ConstIndirectAccessor ConstElemCacheIndirectAccessor;

class Output;

#ifndef ACC
struct neighbor {
    int idx;
    double dist2;
    neighbor() : idx(-1), dist2(0.0) {}
    neighbor(int e, double a) : idx(e), dist2(a) {}
};
#endif

// Update markers in surface elements
struct MarkerUpdate {
    int m;
    int src_elem;
    int dst_elem;
    int inc; // 1=move, 0=remove
    MarkerUpdate() : m(-1), src_elem(-1), dst_elem(-1), inc(0) {}
    MarkerUpdate(int marker, int src, int dst, int increment)
        : m(marker), src_elem(src), dst_elem(dst), inc(increment) {}
};

// Define the struct to store marker data
struct AppendMarkerData {
    double eta[NODES_PER_ELEM];
    int elem;
    int mattype;
    double time;
    double depth;
    double distance;
    double slope;
    int genesis;
    AppendMarkerData() : elem(-1), mattype(-1), time(0.0), depth(0.0), distance(0.0), slope(0.0), genesis(0) {
        for (int i = 0; i < NODES_PER_ELEM; i++) eta[i] = 0.0;
    }
    AppendMarkerData(ConstShapefnAccessor e, int el, int mt, double t, double d, double dis, double s, int g)
        : elem(el), mattype(mt), time(t), depth(d), distance(dis), slope(s), genesis(g) {
        for (int i = 0; i < NODES_PER_ELEM; i++) eta[i] = e[i];
    }
    AppendMarkerData(ConstElemCacheAccessor e, int el, int mt, double t, double d, double dis, double s, int g)
        : elem(el), mattype(mt), time(t), depth(d), distance(dis), slope(s), genesis(g) {
        for (int i = 0; i < NODES_PER_ELEM; i++) eta[i] = e[i];
    }
};

// Update markers in surface elements
struct ElemMarkerInfo {
    int nmarker;
    double coord[NDIMS];
    double value;
    ElemMarkerInfo() : nmarker(0), value(0.0) {
        for (int i = 0; i < NDIMS; i++) coord[i] = 0.0;
    }
    ElemMarkerInfo(const double c[NDIMS], int n, double v, int dst)
        : nmarker(n), value(v) {
        for (int i = 0; i < NDIMS; i++) coord[i] = c[i];
        }
};

typedef std::vector<neighbor> neighbor_vec;
typedef std::vector<MarkerUpdate> MU_vec;
typedef std::vector<AppendMarkerData> AMD_vec;
typedef std::vector<ElemMarkerInfo> EMI_vec;

//
// Structures for input parameters
//
struct Sim {
    double max_time_in_yr;
    double output_time_interval_in_yr;
    int max_steps;
    int output_step_interval;
    int checkpoint_frame_interval;
    int restarting_from_frame;
    int hdf5_compression_level;
    int info_display_step_interval;
    bool is_outputting_averaged_fields;
    bool is_restarting;
    bool has_initial_checkpoint;
    bool has_output_during_remeshing;
    bool has_marker_output;
    int earthquake_output_step_interval;
    bool seismic_moment_calculate_output;
    double earthquake_start_factor;
    double earthquake_end_factor;

    std::string modelname;
    std::string restarting_from_modelname;
};

struct Mesh {
    int meshing_option;
    int meshing_elem_shape;
    int meshing_verbosity;
    bool meshing_sediment;
    int tetgen_optlevel;
    int quality_check_step_interval;

    double xlength, ylength, zlength;
    double resolution;
    double smallest_size;
    double largest_size;
    double sediment_size;
    double min_angle;  // for 2D only
    double min_tet_angle, max_ratio; // for 3D only
    double min_quality;
    double max_boundary_distortion;

    double_pair refined_zonex, refined_zoney, refined_zonez;
    std::string poly_filename;
    std::string exo_filename;

    bool is_discarding_internal_segments;
    int remeshing_option;

    // Parameters for mesh optimizer MMG
    int mmg_debug;
    int mmg_verbose;
    double mmg_hmax_factor;
    double mmg_hmin_factor;
    double mmg_hausd_factor;
    double mmg_init_coarsening_factor;
    bool use_mmg_init;
};

struct Control {
    double gravity;
    double characteristic_speed;
    double inertial_scaling;
    double dt_fraction;
    double fixed_dt;
    double damping_factor;
    int damping_option;
    int ref_pressure_option;
//    bool surface_pressure_correction;
    bool is_using_mixed_stress;
    double mixed_stress_reference_viscosity;

    int surface_process_option;
    double surface_diffusivity;
    double surf_diff_ratio_terrig;
    double surf_diff_ratio_marine;
    double surf_depo_universal;
    double surf_base_level;
    std::string surface_process_gospl_config_file;
    int    gospl_coupling_frequency;  // Run GoSPL every N steps (default: 1)
    double gospl_coupling_interval_in_yr;   // Run GoSPL every T years (default: 1000)
    std::string gospl_coupling_mode;  // "steps" (default) or "time"
    bool gospl_velocity_coupling;  // Send all 3 DES velocity components to GoSPL (default: true)
    double gospl_mesh_resolution;  // GoSPL mesh spacing (default: -1 = auto)
    double gospl_mesh_perturbation;  // Grid perturbation factor 0-1 (default: 0, 0.3=moderate)
    double gospl_mesh_padding;       // Fractional padding beyond DES domain on each side (default: 0.1)
    double terrig_sediment_volume;
    double terrig_sediment_area;
    double terrig_sediment_diffusivity;
    double terrig_depth_coefficient;
    bool is_reporting_terrigenous_info;
    double hemipelagic_sedimentation_rate;
    double hemipelagic_width;
    double hemipelagic_max_depth;
    double pelagic_sedimentation_rate;
    double pelagic_increasing_width;

    bool is_quasi_static;
    bool has_thermal_diffusion;
    bool has_hydraulic_diffusion;

    bool has_hydration_processes;
    double hydration_migration_speed;

    bool has_PT;
    mutable bool PT_jump;
    int PT_max_iter;
    double PT_relative_tolerance;

    bool has_moving_mesh;
    bool use_global_velocity_scaling;

};

struct BC {
    double surface_temperature;
    double mantle_temperature;

    double winkler_delta_rho;
    bool has_winkler_foundation;

    double elastic_foundation_constant;
    bool has_elastic_foundation;

    bool has_water_loading;

    int vbc_x0;
    int vbc_x1;
    int vbc_y0;
    int vbc_y1;
    int vbc_z0;
    int vbc_z1;
    int vbc_n0;
    int vbc_n1;
    int vbc_n2;
    int vbc_n3;

    double vbc_val_x0;
    double vbc_val_x1;
    double vbc_val_y0;
    double vbc_val_y1;
    double vbc_val_x0_l;  // lateral shear value for x0 face (used with bc type 6)
    double vbc_val_x1_l;
    double vbc_val_y0_l;
    double vbc_val_y1_l;
    double vbc_val_z0;
    double vbc_val_z1;
    double vbc_val_n0;
    double vbc_val_n1;
    double vbc_val_n2;
    double vbc_val_n3;

    double bottom_shear_zone_thickness;

    double vbc_val_division_x0_min;
    double vbc_val_division_x0_max;
    double vbc_val_division_x1_min;
    double vbc_val_division_x1_max;

    double vbc_val_x0_ratio0;
    double vbc_val_x0_ratio1;
    double vbc_val_x0_ratio2;
    double vbc_val_x0_ratio3;
    double vbc_val_x1_ratio0;
    double vbc_val_x1_ratio1;
    double vbc_val_x1_ratio2;
    double vbc_val_x1_ratio3;

    int num_vbc_period_x0;
    int num_vbc_period_x1;

    double_vec vbc_period_x0_time_in_yr;
    double_vec vbc_period_x1_time_in_yr;

    double_vec vbc_period_x0_ratio;
    double_vec vbc_period_x1_ratio;

    double vbc_val_z1_loading_period;

    // General stress (Neumann) bcs 
    int stress_bc_x0;
    int stress_bc_x1;
    int stress_bc_y0;
    int stress_bc_y1;
    int stress_bc_z0;
    int stress_bc_z1;

    // hyrdaulic bouncdary
    int hbc_x0;
    int hbc_x1;
    int hbc_y0;
    int hbc_y1;
    int hbc_z0;
    int hbc_z1;

    double stress_val_x0;
    double stress_val_x1;
    double stress_val_y0;
    double stress_val_y1;
    double stress_val_z0;
    double stress_val_z1;
};

struct IC {
    int mattype_option;
    int num_mattype_layers;
    int_vec layer_mattypes;
    double_vec mattype_layer_depths;

    int weakzone_option;
    bool is_restarting_weakzone;
    double weakzone_plstrain;
    double weakzone_azimuth;
    double weakzone_inclination;
    double weakzone_halfwidth;
    double weakzone_y_min;
    double weakzone_y_max;
    double weakzone_depth_min;
    double weakzone_depth_max;
    double weakzone_xcenter;
    double weakzone_ycenter;
    double weakzone_zcenter;
    double weakzone_xsemi_axis;
    double weakzone_ysemi_axis;
    double weakzone_zsemi_axis;
    double weakzone_standard_deviation;
    double weakzone_gaussian_amplitude; // x-shift amplitude for weakzone_option=4 (Gaussian_planar_zone)

    // multi-segment planar weak zone (weakzone_option == 5)
    int weakzone_num_segments;
    double_vec weakzone_segments_xcenter;
    double_vec weakzone_segments_ycenter;
    double_vec weakzone_segments_zcenter;
    double_vec weakzone_segments_azimuth;
    double_vec weakzone_segments_inclination;
    double_vec weakzone_segments_halfwidth;
    double_vec weakzone_segments_x_min;
    double_vec weakzone_segments_x_max;
    double_vec weakzone_segments_depth_min;
    double_vec weakzone_segments_depth_max;
    double_vec weakzone_segments_y_min;
    double_vec weakzone_segments_y_max;

    int temperature_option;
    std::string Temp_filename;
    std::string Nodes_filename;
    std::string Connectivity_filename;
    double oceanic_plate_age_in_yr;
    double continental_plate_age_in_yr;
    double radiogenic_crustal_thickness;
    double radiogenic_folding_depth;
    double radiogenic_heating_of_crust;
    double lithospheric_thickness;
    double radiogenic_heat_dome_center_x;
    double radiogenic_heat_dome_center_y;
    double surface_heat_flux;
    double radiogenic_heat_dome_amplitude;
    double radiogenic_heat_dome_width;
    double radiogenic_heat_dome_width_y;
    double radiogenic_heat_dome_azimuth;
    int nhlayer;
    double_vec radiogenic_heat_boundry;
    int_vec radiogenic_heat_mat_in_layer;

    double isostasy_adjustment_time_in_yr;

    double excess_pore_pressure;
    bool has_body_force_adjustment;
};

struct Mat {
    int rheol_type;
    int phase_change_option;
    int nmat;
    int mattype_ref;
    int mattype_mantle;
    int mattype_depleted_mantle;
    int mattype_partial_melting_mantle;
    int mattype_crust;
    int mattype_sed;
    int mattype_oceanic_crust;
    int mattype_mor_extrusion;
    int mattype_asthenosphere;
    double convert_rate_oceanic_crust;

    bool is_plane_strain;
    double visc_min;
    double visc_max;
    double tension_max;
    double therm_diff_max;

    double_vec rho0;
    double_vec alpha;

    double_vec bulk_modulus;
    double_vec shear_modulus;

    double_vec visc_exponent;
    double_vec visc_coefficient;
    double_vec visc_activation_energy;
    double_vec visc_activation_volume;

    double_vec heat_capacity;
    double_vec therm_cond;
    double_vec radiogenic_heat_prod;

    // plastic parameters
    double_vec pls0, pls1;
    double_vec cohesion0, cohesion1;
    double_vec friction_angle0, friction_angle1;
    double_vec dilation_angle0, dilation_angle1;

    // hydraulic parameters
    double_vec porosity;
    double_vec hydraulic_perm;
    double_vec fluid_rho0;  // pore fluid density
    double_vec fluid_alpha; // pore fluid thermal expansivity
    double_vec fluid_bulk_modulus;  // pore fluid bulk modulus
    double_vec fluid_visc;  // pore fluid dynamic viscosity
    double_vec biot_coeff;  // Biot-Willis coefficient
    double_vec bulk_modulus_s;  // bulk modulus of solid grain (mineral)
  
    // rate-and-state friction parameters
    double_vec direct_a;
    double_vec evolution_b;
    double_vec characteristic_velocity;
    double_vec characteristic_distance;
    int state_var_model;
    // double_vec static_friction_coefficient;

};

struct Time {
    int64_t remesh_time;
    int64_t output_time;
    int64_t start_time;
};

struct Markers {
    int init_marker_option;
    int markers_per_element;
    int min_num_markers_in_element;
    int replenishment_option;
    uint random_seed;
    double init_marker_spacing;
};

struct Debug {
    bool dt;
//    bool has_two_layers_for;
};

enum MonitorRebindMode {
    monitor_rebind_initial_coord = 0,
    monitor_rebind_pre_remesh_coord = 1
};

struct Monitor {
    bool enabled;
    int step_interval;
    int num_points;

    std::string points_unit;
    double points_scale_to_m;
    int remesh_rebind_mode;

    double_vec points_x;
    double_vec points_y;
    double_vec points_z;

    std::string output_prefix;
    bool write_header;

    bool output_coord;
    bool output_velocity;
    bool output_force;
    bool output_temperature;
    bool output_pore_pressure;
    bool output_bcflag;

    bool output_stress;
    bool output_strain;
    bool output_strain_rate;
    bool output_plastic_strain;
    bool output_plastic_strain_rate;

    bool output_radiogenic_source;
    bool output_density;
    bool output_mesh_quality;
    bool output_viscosity;
    bool output_material;

    bool output_dynamic_friction;
    bool output_state_variable;
};

struct PointCloud {
    const array_t &data;

    PointCloud(const array_t &_data)
        : data(_data) {}

    inline size_t kdtree_get_point_count() const { return data.size(); }

    inline double kdtree_get_pt(const size_t idx, const size_t d) const {
        return data[idx][d];
    }

    template <class BBOX>
    bool kdtree_get_bbox(BBOX&) const { return false; }
};

using NANOKDTree = nanoflann::KDTreeSingleIndexAdaptor<nanoflann::L2_Simple_Adaptor<double, PointCloud>, PointCloud, NDIMS>;

struct Param {
    Sim sim;
    Monitor monitor;
    Mesh mesh;
    Control control;
    BC bc;
    IC ic;
    Mat mat;
    Markers markers;
    Debug debug;
};

//
// Structures for surface processes
//
struct SurfaceInfo {

//    const double sec_year = 31556925.2;

    int efacet_top;
    int ntop;
    int etop;

    double base_level;
    double surf_diff;
    double diff_ratio_terrig;
    double diff_ratio_marine;
    double terrig_diffusivity;
    double terrig_dpeth_coeff;
    double depo_universal;
    double ero_rate;
    double max_surf_vel;

    double_vec *dh;
    double_vec *dh_oc;
    double_vec *src_locs;

    double_vec *total_dx;
    double_vec *total_slope;

    int_vec *top_nodes;
    int_vec *top_facet_elems;

    int_vec *landform_map;
    double_vec *drainage;

    double_vec *dhacc;
    // variable allocate by remesh
    double_vec *edvacc_surf;
    int_vec2D *node_and_elems;
    segment_t *elem_and_nodes;

    int_map arctop_facet_elems;
    int_map arctop_nodes;


    int ntops;
    int nbots;
    int_vec tops;
    int_vec coasts;
    int_vec bots;

    double_vec top_elev;
    double_vec bot_elev;

};

//
// Structures for model variables
//
class MatProps;
class MarkerSet;
struct Variables {
    double time;
    double dt;
    double dt_PT;
    double l2_residual;
    int steps;
    int nremesh;
    int noutput;
    Time func_time;
    Output *output;

    int nnode;
    int nelem;
    int nseg;
    int ncontact;
    int nx, ny, nz, ncell;

    double max_vbc_val;
    double max_global_vel_mag;
    double global_dt_min;
    double compensation_pressure;
    double bottom_temperature;

// #ifdef ATS
//     double vmax, hmin, dt_elastic, dt_min, vmax_shear_zone, CL_min;
// #endif    
//     double a0, b0, a1, b1, max_pls, number_plf, seismic_eff, S_E, K_E;

//     double_vec *MAX_shear, *CL, *dl_min, *maxv, *strain_energy, *kinetic_energy, *MAX_shear_0;
// #ifdef RS
//     double_vec *state1, *slip_velocity, *friction_coefficient, *RS_shear, *Failure_mode, avg_shear_stress, avg_vm, slip_area;
// #endif

    // These 5 arrays are allocated by external library
    array_t *coord;
    conn_t *connectivity;
    conn_t *connectivity_surface;
    segment_t *segment;
    segflag_t *segflag;
    regattr_t *regattr;
    regular_t *cell;


    uint_vec *bcflag;
    int_vec *bnodes[nbdrytypes];
    std::vector< std::pair<int,int> > *bfacets[nbdrytypes];
    array_t *bnormals;
    int vbc_types[nbdrytypes];
    int hbc_types[nbdrytypes_hydro];
    int stress_bc_types[nbdrytypes_hydro];
    double vbc_values[nbdrytypes];
    double stress_bc_values[nbdrytypes];
    double vbc_val_z1_loading_period;

    std::map<std::pair<int,int>, double_vec> edge_vectors;
    double_vec edge_vec;
    int_vec edge_vec_idx;
    double_vec vbc_vertical_div_x0;
    double_vec vbc_vertical_div_x1;
    double_vec vbc_vertical_ratio_x0;
    double_vec vbc_vertical_ratio_x1;

    int_vec *top_elems;
    int ntop_elems;
    int_vec2D *markers_in_elem;
    int_vec2D *hydrous_markers_in_elem;

    int_vec2D *support;
    int_vec *support_arr;
    int_vec *support_idx;
    conn_t *neighbor; // neighboring elements for each element
    int_pair_vec *contact; // contact elements for each element
    double_vec *ctmp; // temporary array for contact elements

    double_vec *volume, *volume_old, *volume_n;
    double_vec *mass, *tmass;
    double_vec *hmass;
    double_vec *ymass; // Young's modulus for nodes
    double_vec *edvoldt;
    double_vec *temperature, *plstrain, *delta_plstrain;
    double_vec *stressyy, *dpressure, *viscosity;
    double_vec *old_mean_stress;
    double_vec *ntmp;
    double_vec *init_elem_size_n;   // frozen initial nodal element size for MMG metric
    double_vec *radiogenic_source;

    // For hyraulic proceses
    double_vec *fmass; // pore water mass
    double_vec *ppressure; // pore pressure
    double_vec *dppressure; // delta pore pressure
    double_vec *dppressure_zero; // delta pore pressure
    double_vec *fluid_source; // injection and pumping of pore water

    // For rate-and-state friction
    double_vec *dyn_fric_coeff;
    double_vec *state_variable;
    
    // For surface processes
    SurfaceInfo surfinfo;
    int_vec melt_markers;

#ifdef HAS_GOSPL_CPP_INTERFACE
    // GoSPL driver for landscape evolution modeling
    GoSPLDriver* gospl_driver;
#endif

    array_t *vel, *force, *coord0;
    array_t *force_residual;
    tensor_t *strain_rate, *strain, *stress;
    elem_cache *tmp_result;
    double_vec *etmp;
    int_vec *etmp_int;

    // tensor_t *stress_old;

    MatProps *mat;

    std::vector<MarkerSet*> markersets;
    int hydrous_marker_index;

    int_vec2D *elemmarkers; // for marksersets[0] (mattype markers)
    Array2D<int,1> *hydrous_elemmarkers; // for markersets[hydrous_marker_index] (hydrous markers)

    Variables()
    {
        vbc_vertical_div_x0.resize(4);
        vbc_vertical_div_x1.resize(4);
        vbc_vertical_ratio_x0.resize(4);
        vbc_vertical_ratio_x1.resize(4);
#ifdef HAS_GOSPL_CPP_INTERFACE
        gospl_driver = nullptr;
#endif
    }

    double_vec *log_table; // for log(x) lookup table
    double_vec *tan_table; // for tan(x) lookup table
    double_vec *sin_table; // for sin(x) lookup table

};

#endif