Calculating double derivatives options

Hi, In Paraview, if I want to calculate double derivative/second derivative of a variable, we have got the following options: 1) Take gradient filter twice-but this is associated with loss of accuracy and can produce some CFD mesh artefacts. Can someone please let me know if whether adding ghost cells can remove these artefacts? And how to use ghost cells? 2) Can i write double derivative code/function from theory in the programmable filter? If so, can some one let me know how to do it? Or do we have to export to MATLAB or Python if we want to calculate the double derivative efficiently? Any suggestions please.thanks

Calculating higher-order precision gradients using vmtk’s vtkvmtkUnstructuredGridGradientFilter in Python is relatively straightforward.

For an Anaconda environment, first create a virtual environment and install vmtk using the following commands:

conda create -n vmtk python=3.11
conda activate vmtk
conda install conda-forge::vmtk

Next, you can create a Hessian matrix by applying vtkvmtkUnstructuredGridGradientFilter twice to a variable with a Python program like the one below. However, this program requires modification as it takes a vtu file format as input. You will likely also need to customize the arguments and other settings.

#!/usr/bin/env python3
"""
VMTK Higher-Order Precision Gradient Calculator

This script calculates first and second-order gradients of a specified field
in a VTK unstructured grid file using VMTK's differential geometry filters.
"""

import argparse
import os
import sys
import vtk
from vmtk.vtkvmtkDifferentialGeometryPython import vtkvmtkUnstructuredGridGradientFilter


def parse_arguments():
    """Parse command line arguments."""
    parser = argparse.ArgumentParser(
        description="This program computes the Hessian matrix using vmtk's higher-order precision gradient calculation filter",
        formatter_class=argparse.ArgumentDefaultsHelpFormatter
    )
    parser.add_argument(
        '--input-file', 
        default='test.vtu', 
        help='Input VTK unstructured grid file'
    )
    parser.add_argument(
        '--output-file', 
        default='out.vtu', 
        help='Output file for gradient results'
    )
    parser.add_argument(
        '--input-array-name', 
        default='DistanceToCenter', 
        help='Scalar field to calculate gradients from'
    )
    parser.add_argument(
        '--quadrature-order', 
        type=int, 
        default=3, 
        help='Quadrature order for gradient calculations'
    )
    return parser.parse_args()


def calculate_gradients(input_file, array_name, quadrature_order):
    """
    Calculate first and second order gradients.
    
    Args:
        input_file: Path to input VTU file
        array_name: Name of the scalar field for gradient calculation
        quadrature_order: Order of quadrature for numerical integration
        
    Returns:
        VTK data object with calculated gradients
    """
    # Check if input file exists
    if not os.path.exists(input_file):
        raise FileNotFoundError(f"Input file not found: {input_file}")
    
    # Read the input file
    reader = vtk.vtkXMLUnstructuredGridReader()
    reader.SetFileName(input_file)
    reader.Update()
    
    # Calculate first-order gradient
    gradient_filter1 = vtkvmtkUnstructuredGridGradientFilter()
    gradient_filter1.SetInputConnection(reader.GetOutputPort())
    gradient_filter1.SetInputArrayName(array_name)
    gradient_filter1.SetGradientArrayName(f'{array_name}_grad1')
    gradient_filter1.SetQuadratureOrder(quadrature_order)
    
    # Calculate second-order gradient
    gradient_filter2 = vtkvmtkUnstructuredGridGradientFilter()
    gradient_filter2.SetInputConnection(gradient_filter1.GetOutputPort())
    gradient_filter2.SetInputArrayName(f'{array_name}_grad1')
    gradient_filter2.SetGradientArrayName(f'{array_name}_grad2')
    gradient_filter2.SetQuadratureOrder(quadrature_order)
    gradient_filter2.Update()
    
    return gradient_filter2.GetOutput()


def save_results(data, output_file):
    """
    Save the VTK data to an output file.
    
    Args:
        data: VTK data object to write
        output_file: Path to output file
    """
    # Ensure the output directory exists
    output_dir = os.path.dirname(output_file)
    if output_dir and not os.path.exists(output_dir):
        os.makedirs(output_dir)
        
    # Write the output file
    writer = vtk.vtkXMLUnstructuredGridWriter()
    writer.SetFileName(output_file)
    writer.SetInputData(data)
    writer.Write()
    
    print(f"Results successfully written to {output_file}")


def main():
    """Main execution function."""
    try:
        # Parse command line arguments
        args = parse_arguments()
        
        # Calculate gradients
        result_data = calculate_gradients(
            args.input_file,
            args.input_array_name,
            args.quadrature_order
        )
        
        # Save results
        save_results(result_data, args.output_file)
        
    except Exception as e:
        print(f"Error: {e}", file=sys.stderr)
        return 1
        
    return 0


if __name__ == "__main__":
    sys.exit(main())

vmtk_test.zip (88.7 KB)