Morphogenesis is a highly complex multi-scale phenomenon governed by the interplay between intrinsic biological processes and external fluid dynamics, involving the dynamic remodeling of biological structures under hydrodynamic stress. Modeling the topological evolution of such developing substrates presents significant numerical challenges. Conventionally, phase-field models like the Cahn-Hilliard equation are employed for these diffuse interface problems; however, their fourth-order spatial derivatives impose severe time-step restrictions and require wide computational stencils, creating a significant bottleneck for long-time simulations. To overcome this limitation, we present an integrated computational framework designed to simulate active morphogenesis using a conservative Allen-Cahn approach. By reducing the governing PDE to a second-order reaction-diffusion system, we significantly enhance computational efficiency while rigorously maintaining mass conservation and consistent interfacial forces. Adopting a "one-fluid" continuous approach, the substrate is modeled as a developing viscoelastic continuum coupled with the incompressible Navier-Stokes equations. The framework accounts for the complex rheology of the biological material using an Oldroyd-B constitutive model, while growth is driven by the advection-diffusion-reaction of a nutrient scalar, creating a fully coupled metabolic feedback loop. The numerical solution employs finite-difference schemes on Cartesian grids and a projection method, coupled with high-order Implicit-Explicit (IMEX) time integration to efficiently resolve stiff terms and high viscosity contrasts. Particular attention is devoted to the discretization of the coupled system, with emphasis on mass conservation, boundedness of the phase field, and the consistent transfer of stresses between the phase-field and momentum equations. Finally, we demonstrate the method's capabilities through examples of substrate morphogenesis within complex flow regimes, including surface-attached biofilm formation. The framework is benchmarked against published solutions to quantify deviations in thermodynamic consistency, highlighting the favorable trade-off between algorithmic efficiency and physical fidelity in dynamic fluid-structure interaction scenarios.