I doubted this elbo function. Usually, in the examples, we have the loss function implemented in the following way:
def loss_function(recon_x, x, mu, logvar):
    BCE = F.binary_cross_entropy(recon_x, x.view(-1, 784), reduction="sum")

    # see Appendix B from VAE paper:
    #  Kingma and Welling. Auto-Encoding Variational Bayes. ICLR, 2014
    # https://arxiv.org/abs/1312.6114
    # 0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())

    return BCE + KLD, BCE, KLD

I wanted to reproduce this using the torch.nn.functional.kl_div function, which is similar to what OP implemented in elbo:
def loss_function(recon_x, x, mu, logvar):
    BCE = F.binary_cross_entropy(recon_x, x.view(-1, 784), reduction="sum")

    # see Appendix B from VAE paper:
    # Kingma and Welling. Auto-Encoding Variational Bayes. ICLR, 2014
    # https://arxiv.org/abs/1312.6114
    # 0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    
    KLD = torch.nn.functional.kl_div(
        recon_x, 
        x.view(-1, 784), 
        reduction='batchmean', 
        log_target=False
    )

    return BCE + KLD, BCE, KLD

But, I wonder — Is this correct? I ask that because the KLD should measure the divergence between the mean of the posterior (latent distribution) with respect to the mean of the prior (normal distribution). One can see that this does not work because reconstructions are sharp, but the generative aspect of the VAE is completely lost compared to the loss when you use
KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp()).
I think that kl_div should be receiving mu that comes from the latent distribution, and for the target, it should be receiving the eps from the reparameterize function?
def reparameterize(self, mu, logvar):
    
    std = torch.exp(0.5 * logvar)
    eps = torch.randn_like(std)

    return mu + eps * std, eps

PS. I probably should have opened a new post for this, but it made sense here in my mind.