Skip to content

Custom AutoregressiveCovariance operator #5018

@nemaniravi7

Description

@nemaniravi7

Anisotropic Background covariance operator
Pass different length scales, with different tensor shapes with 2D or 3D space to model the covariances.
I can pass the diffusion form but I can't pass the custom weight through firedrake's inbuilt function.

mesh = RectangleMesh(20, 10, 2.0, 1.0)
V = FunctionSpace(mesh, "CG", 1)
 x, y = SpatialCoordinate(mesh)
    
 u0_expr = exp(-50 * ((x - 1.0)**2 + (y - 0.5)**2))
innovation = Function(V).interpolate(u0_expr)

# Use Firedrake's kappa_m to get exact diffusion coefficients
        kappa_x = kappa_m(Constant(Lx), m_steps)
        kappa_y = kappa_m(Constant(Ly), m_steps)
            
# Build the Anisotropic Tensor using the computed kappas
 theta = Constant(np.radians(angle_deg))
 c, s = cos(theta), sin(theta)
 R_mat = as_matrix([[c, -s], [s, c]])
 D_mat = as_matrix([[kappa_x, Constant(0.0)], [Constant(0.0), kappa_y]])
 K_tensor = dot(R_mat, dot(D_mat, R_mat.T))
        
 u = TrialFunction(V)
 v = TestFunction(V)
        
 # Explicitly add the Mass Term to prevent DIVERGED_LINEAR_SOLVE
 custom_diffusion_form = inner(u, v) * dx + inner(dot(K_tensor, grad(u)), grad(v)) * dx
        
 # Instantiate the official operator
 B_operator = AutoregressiveCovariance(V, Lx, sigma=1.0, m=m_steps, form=custom_diffusion_form)
        
 # Pannekoucke & Massart (2008) Anisotropic Gamma Normalization #https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.288)
 d = 2 
 gamma_ratio = math.gamma(m_steps / 2) / math.gamma((m_steps / 2) - (d / 2))
 spatial_scale = (4.0 * math.pi) * (Lx * Ly) 
 lam = spatial_scale * gamma_ratio
        
 B_operator._weight = Constant(lam) # custom weight to bypass Firedrake's internal _weight attribution

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions