LLNL / conduit

Simplified Data Exchange for HPC Simulations
https://software.llnl.gov/conduit/
Other
207 stars 63 forks source link

Add Gyre as example dataset #1257

Open nicolemarsaglia opened 4 months ago

nicolemarsaglia commented 4 months ago

gyre ref: https://shaddenlab.berkeley.edu/uploads/LCS-tutorial/examples.html

void
tutorial_gyre_example(float64 time_value, Node &mesh)
{
    mesh.reset();
    int xy_dims = 40;
    int z_dims = 2;

    conduit::blueprint::mesh::examples::braid("hexs",
                                             xy_dims,
                                             xy_dims,
                                             z_dims,
                                             mesh);

    mesh["state/time"] = time_value;
    Node &field = mesh["fields/gyre"];
    field["association"] = "vertex";
    field["topology"] = "mesh";
    field["values"].set(DataType::float64(xy_dims*xy_dims*z_dims));

    Node &vec_field = mesh["fields/gyre_vec"];
    vec_field["association"] = "vertex";
    vec_field["topology"] = "mesh";
    vec_field["values/u"].set(DataType::float64(xy_dims*xy_dims*z_dims));
    vec_field["values/v"].set(DataType::float64(xy_dims*xy_dims*z_dims));
    vec_field["values/w"].set(DataType::float64(xy_dims*xy_dims*z_dims));

    float64 *values_ptr = field["values"].value();
    float64 *u_values_ptr = vec_field["values/u"].value();
    float64 *v_values_ptr = vec_field["values/v"].value();
    float64 *w_values_ptr = vec_field["values/w"].value();

    float64 e = 0.25;
    float64 A = 0.1;
    float64 w = (2.0 * PI_VALUE) / 10.0;
    float64 a_t = e * sin(w * time_value);
    float64 b_t = 1.0 - 2 * e * sin(w * time_value);
    // print("e: " + str(e) + " A " + str(A) + " w " + str(w) + " a_t " + str(a_t) + " b_t " + str(b_t))
    // print(b_t)
    // print(w)
    int idx = 0;
    for (int z=0; z < z_dims; z++)
    {
        for (int y=0; y < xy_dims; y++)
        {
            // scale y to 0-1
            float64 y_n = float64(y)/float64(xy_dims);
            float64 y_t = sin(PI_VALUE * y_n);
            for (int x=0; x < xy_dims; x++)
            {
                // scale x to 0-1
                float64 x_f = float(x)/ (float(xy_dims) * .5);
                float64 f_t = a_t * x_f * x_f + b_t * x_f;
                // print(f_t)
                float64 value = A * sin(PI_VALUE * f_t) * y_t;
                float64 u = -PI_VALUE * A * sin(PI_VALUE * f_t) * cos(PI_VALUE * y_n);
                float64 df_dx = 2.0 * a_t + b_t;
                // print("df_dx " + str(df_dx))
                float64 v = PI_VALUE * A * cos(PI_VALUE * f_t) * sin(PI_VALUE * y_n) * df_dx;
                values_ptr[idx] = sqrt(u * u + v * v);
                u_values_ptr[idx] = u;
                v_values_ptr[idx] = v;
                w_values_ptr[idx] = 0;
                // values[idx] = u * u + v * v
                // values[idx] = value
                // print("u " + str(u) + " v " + str(v) + " mag " + str(math.sqrt(u * u + v * v)))
                idx++;
            }
        }
    }