I would claim a valid runtime error raised by an invalid mathematical operation is way less painful than an unexpected behavior by the framework, which would change the behavior of the model and could potentially diverge your training.
Debugging such divergence issues can take many more weeks than to fix a shape issue.
As described before, I would see the workaround using the nn.Identity module as a silent numerical error, since the output would not be normalized at all and could have a different range and thus cause divergence.
Here is a small example:
bn = nn.BatchNorm1d(3)

x = torch.randn(2, 3, 1) * 1000.
out = bn(x)
print(out.min(), out.max(), out.mean())
# > tensor(-1., grad_fn=<MinBackward1>) tensor(1., grad_fn=<MaxBackward1>) tensor(0., grad_fn=<MeanBackward0>)

# workaround with Identity 
bn = nn.Identity()
out = bn(x)
print(out.min(), out.max(), out.mean())
# > tensor(-1368.4515) tensor(1302.6108) tensor(75.8487)

x is the incoming activation with a large value range (randn * 1000.). In the working case the batchnorm layer would normalize it and the following layer (e.g. a linear layer) would expect to see these normalized values. If you are now skipping the norm the following layer would get the unnormalized values with an entirely different range.
You could see similar effects by using a pretrained model (e.g. use torchvision.models.resnet18) and by comparing the achieved accuracy using normalized input images (expected) vs. unnormalized inputs with values in the range [0, 255]. In the latter case the model would not perform well even though it’s a proper pretrained model.